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SEISMIC  DIFFRACTION  BY  NON-UNIFORl  CRUSTAL  STRUCTURES 


by  Julius  Kane 


abstract 

The  final  report  on  contract  AF  19(628)  -  319  describes  the 
research  contr  Voutions  of  a  group  studying  theoretical  problems  of  sela> 
mology  at  the  University  of  Rhode  Island.  These  are  described  in  simp¬ 
lified  terms  and  formulations  which  are  mathematically  tractable  and 
still  relevent  to  realistic  crustal  structures.  Reprints  of  papers 
describing  the  methods  are  included  as  appendices  in  the  report. 
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SUMMARY  OF  WORK  DOFE  U.NDSR  SUBJECT  COIiTRACT 

The  mathematical  problems  of  theoretical  seismology  are 
especially  acute  for  two  distinct  wave  equations  need  be  satisfied. 
Furthermore,  at  refracting  interfaces,  the  continuity  relations 
usually  demand  that  complex  transcendental  relationships  be  satisfied. 

For  these  reasons,  the  geometric  configurations  for  which  seismic 
problems  are  tractable  are  quite  few.  Separable  geometries  include 
parallel  stacks  of  plane  layers,  and  concentric  spherical  shells  of 
differing  elastic  media.  The  separability  of  the  wave  equations  in 
the  latter  case  is  particularly  fortunate  because  it  describes  a 
spherical  geometry  which  duplicates  the  actual  situation  quite 
faithfully.  However,  geophysical  waveguides  encountered  in  practice 
contain  many  irregularities  whose  variation  does  not  fall  within  the 
scope  of  separable  analysis.  Such  configurations  include  transition 
regions  at  the  continental  margins,  and  tapering  crustal  layers.  In 
addition,  unexplained,  and  possibly  important  diffraction  effects 
arise  from  wall  irregularities;  that  is,  the  interfaces  bounding 
geophysical  waveguides  are  not  smooth,  but  rough  in  nature.  In  the 
past,  such  effects  could  have  been  neglected,  but  with  the  desire  of 
the  seismologist  to  continually  resolve  finer  details  of  the  earth's 
crustal  structure,  such  effects  can  no  longer  be  ignored.  High 
resolution  means  that  wc  need  analyse  and  understand  the  high  frequency 
portions  of  seismic  records.  Interpreting  such  data  in  terms  of  abstract 
idealisations  can  be  misleading  and  deceive  the  geologist  in  his 
understanding  pf  the  subterranean  characteristics  of  the  world. 

J* 

Questions  of  this  nature  have  prompted  our  investigations  under  the 
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of  the  above  problem  in  that  ve  studied  a  layered  elastic  surface 
vhich  is  bent  through  an  angle.  In  this  geometry,  the  layering 
recalls  that  an  infinite  number  of  inodes  contribute  to  the  field'. 

Of  course,  for  any  given  frequency,  only  a  fin_tc  number  of  these 
are  propagating  modes.  To  make  the  problem  a  little  more  interesting, 
vc  considered  the  case  vhen  tvo  modes  can  propagate.  The  methods  ve 
used  for  this  problem  were  not  complex  variable-  techniques,  but  rather 
tools  from  the  theory  of  functional  analysis.  That  is,  we  matched 
fields  on  tve  sides  of  some  conveniently  chosen  "aperture  plane," 

In  terms  of  a  suitable  metric,  a  least  squares  fitting  is  automatically 
satisfied.  Mathematically  speaking,  this  is  an  interesting  method 
because  it  is  applicable  to  a  vide  variety  of  otherwise  intractable 
problens.  Included  within  the  material  of  Appendix  B  are  various 
graphs  showing  the  reflection,  transmission,  and  conversion  coefficients 
as  functions  of  wedge  angle  and  frequency.  One  conclusion  not 
previously  established  is  that  seismic  energy  has  a  great  tendency  to 
continue  propagating  in  the  same  direction,  preferring  to  be  converted 
into  a  different  forward  mode  rather  than  a  reflected  one  of  the  same 
type  • 

During  the  arse  of  our  studies,  it  became  increasingly  evident 
that  solution  of  g  eralized  problems  depended  upon  tbe  solution  of 
the  second  order  linear  differential  equation  with  coefficients  xhat 
might  be  substantially  different  from  the  standard  cases  treated  in 
the  literature.  For  this  reason  we  considered  the  relevent  mathematical 
problem  from  a  computational  point  of  view,  one  which  would  furnish 


-  6  - 


^bcrs  and  graphs  in  a  atraight-forvard  manner  vlthout  being  M 
by  ques lions  of  turning  points  and  transition  regions.  The  raault.  of 
this  are  presented  in  Appendix  C.  I^n^of  tho  Second-order 

-  m.-Orcitial  So.ua felon  vhich  was  published  in  tho  Journal  of 

Mat'-. era  ti.cai  Fnysic  s . 

An  immediate  application  of  the  results  mentioned  in  the 
paragraph  above  was  a  vay  of  estimating  the  corruption  of  tine/phase 
data  by  random  inhdmo5enelties  along  the  propagation  path.  Assuming 
only  that  the  mean  square  variation  of  the  velocity  profile  is  kroon, 

„e  were  able  to  find  a  correction  formula  for  dispersion  characteristics 
Cc  waves  propagating  along  inh«Wgenooua  paths.  Tb.  formula  given  is 
a  very  simple  one  and  Is  given  expllcity  as  equation  (27)  of  App-ndi*  D. 

I„  Appendix  E  *e  found  a  method  to  use -dispersion  relations 
between  interacts  vaveguides  as  a  novel  type  of  strain  gunge.  In 
this  ease,  dispersion  of  optical  surface  oaves  la  used  to  measure^  ^ 
lengths  or  strains  to  high  accuracy,  of  the  order  of  one  part  In  lO  . 

In  Appendix  F,  oe  present  a  neo  way  of  solving  diffraction 
problems.  Customarily,  the  technique  has  been  to  distort  the  geometry 
0r  crustal  configurations  and  consider  teem  to  be  of  idealined,  but 
highly  distorted  types:  step  function  changes  in  elevation  and  other 
abrupt  transitions.  In  actual  crustal  configurations,  staircase 
elevations  arc  not  encountered,  vaveguides  are  joined  smoothly  by 
60M  transition  region.  Bather  than  pose  a  boundary  value  problem, 
lt  seemed  reasonable  to  ask  the  question:  In  vhat  vays  can  seismic 
energy  be  converted  from  the  modes  of  one  structure  »  another!  Vc 
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expected  that  the  answer  would  not  oe  unique.  After  all,  there  are 
a  nondenumerab-e  infinity  oi  ways  by  which  the  t’-nsition  can  be 
effected.  On  the  other  hand,  we  did  expect  that  the  low  order 
effects  would  remain  stable;  that  is,  relatively  insensitive  to  the 
type  of  transition  chosen.  This,  of  course,  is  the  principle  behind 
assuring  step  function  changes.  However,  in  our  analysis,  we  were 
not  confined  to  step  function  changes,  and  the  type  transition  we 
discussed  could  easily  be  extended  to  cover  the  case  of  more  involved 
geometries.  The  future  of  this  method  seems  particularly  promising, 
especially  in  light  of  the  ease  by  which  the  results  in  Appendix  F 
were  found.  Nothing  more  than  the  solution  of  five  algebraic 
equations  in  five  unknowns  was  required,  and  even  that  required  no 
matrix  algebra  because  a  simple  substitution  allowed  the  equations 
to  be  solved  by  inspection.  After  the  termination  of  the  present 
contract,  subsequent  studies  have  been  initiated.  Based  upon  our 
preliminary  results,  we  can  mention  that  it  now  seems  as  though 
diirraction  by  arbitrary  transitions  between  geophysical  waveguides 
are  easily  handled. 

Finally,  in  Appendix  G  we  describe  the  results  of  a  computer 
study  of  the  array  problem.  Considering  the  large  investment  embodied 
in  such  projects  as  the  LASA  array,  it  seemed  reasonable  to  ask  the 
following  question;  How  severe  might  .the  effects  of  irregular  crustal 
layering  be  on  its  discrimination  characteristics?  Already,  some 
field  estimates  have  been  given  that  indicate  the  theoretical  signal 
to  noise  ratio  falls  below  expectations.  In  view  of  this,  ve 


in  a  leered  mo  rvrtvr.nvcion  e«.cets  oi  a  at&cx  ox  saew  lavers  on 


sci?:.<ic  fu"  ■".  This  is  or."  p-ob  ra  vh  -oi :.  the  hai-.-_ur.ic  problsa 
g  — g  r.o.,  to-.. .  an  a. «  ’nuit.i. .  -1.(1  -  * .  -*-.  i  ,  y  i.'i*  uj,  in.  or.'.uiio-.i 

?ox*  that  reason  vc  worked  dirc-o  .fy  in  the  tine  dona  in,  thus  sparing 
ourselves  or  the  burden  or  wiving  extensive  Fourier  transforms .  In 
Append!.*;  G,  wo  sus:  arise  the  distortion  of  a  plane  pulse  as  it 
propagates  through  a  uniform  dipping  crust.  The  results  show  that 
distortion  across  the  elements  of  the  array  is  to  he  expected,  and, 
o"ing  to  reverberation  effects,  this  distortion  increases  with  time . 
However,  there  are  simple  euqiization  techniques  which  can  restore 
the  signals  to  a  standard  re*  srer.ee.  The  "rubber  sheet  equalisation" 
is  fully  discussed  in  Appendix  C. 

Appendices  H  and  I  susmarlse  sc-  e  of  the  last  work  done  on 
the  contract.  In  the  last  norths ,  ve  followed  up  the  work  done  in 
Appendix  G  to  ai-ow  for  stacks  of  arbitrary  skew  layers.  In  addition, 
ve  initiated  sene  work  or.  aero-frequency  seismology  —  the  static 
problem.  In  Appendix  I  we  sursaauize  a  novel,  approach  to  the 
clastos vatic  problem.  Reprints  of  there  papers  will  ce  forwarded  to 


he  contracting  officers  once  they  arc  uvuijaoie, 


Appendix  J  is  an  abstract  of  a  master's  t.;csi«  ’  u ncc* 

this  contract,  but  r.ot  published.  Mr.  David  Defanti's  res.:- to  extend 
and  elaborate  upon  some  of  the  published  results  discussed  in  the 
urc vices  abstracts. 
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! 

vf  ;d  upon  our  experience  and  published  results,  we  feel  that 
the  techniques  of  analytic  continuation  can  contribute  a  great  deal 
to  solving  the  problems  of  theoretical  seismology,  particularly  those 
dealing  with  diffraction  by  realistic  crustal  structures  *  The  chief 

i 

advantage  of  this  method  i.i  that,  unlike- separation  of  variables,  :it 
is  not  confined  to  standard  geometries.  We  hope  to  be  able  to 
develop  these  methods  further  at  a  future  date. 

Much  work  still  needs  to  be  done  in  time-deperdent  diffraction. 
While  Fourier  analysis  has  often  been  proposed  as  a  means  to  translate 
harmonic  information  inti  time  domain  information,  tnis  is  in  practice 
seldom  done.  The  reasor,  of  course,  is  the  large  amount  of  computer 
time  to  take  the  required  Fourier  transforms.  Furthermore,  there  are 
theoretical  difficulties:  For  example,  the  spectrum  of  actual  pulses 
is  rot  time- invar lent. 

In  the  same  fashion,  the  operation  of  seismic  arrays  -ir.  the 
time  dor.ain  needs  to  be  fvuly  explored.  While  array  theory  has  been 
extensively  developed  for  electromagnetic  antennas,  the  seismic  problem 
is  quilt  different:  Transients  play  n  bigger  role,  and  the  propagation 
taedium  is  anything  hut  homogeneous.  Also,  in  electromagnetic  problems 
reverberations  can  be  ignored,  but  this  ia  not  the  situation  in  the 
seismic  case. 

Any  inv  rse  problem  is  highly  unstable,  and  approaches  suCh  as 
those  in  Appendix  I  should  be  followed  up.  a*.  effect,  what  we  did 
was  to  admit  that  the  complete  characteristics  of  the  source  c  ;ld  not 
be  recovered  with  perfect  fidelity.  But  even  granting  this,  there 
is  still  the  question:  Can  ve  say  anything  about  the  source  with 


STUDENT  AND  FACULTY  CONTRIBUTORS 


pernops  even  rare  important  than  the  vritten  reports  ana  research 
per formed  under  contract  sponsorship  arc  the  contributions  t.  at 
will  be  made  by  student  and  faculty  contributors  whose  further 
research  participation  training  has  been  made  possible  in  part 
by  the  contract's  financial  support.  Personnel,  in  addition  to  the 
principal  investigator  Dr.  Julius  Kane,  include: 


Faculty 


t 

Dr.  John  5per.ee,  Assistant  Professor  of  Electrical  Engineering 


"  _  Suryanaraynn,  Associate  Professor  of  Mathematics 

at  U.R.I. 

Dr.  A.  Mai,  Assistant  Research  Geophysicist  at  U.C.LA. 

Dr.  D.  Sadeh,  Assistant  Research  Geophysicist  at  U.C.L.A. 

Students 

Mr.  David  pence,  graduate  student  at  U.R.I. 

Mr.  Stewart  Davis,  graduate  student  at  U.R.I. 

Mr.  J.  M  score,  graduate  student  at  U.R.I. 

Mr.  Jose  Chirivella,  graduate  student  at  Caltech. 

Mr.  John  Rouse,  student  at  Caltech. 

Mr.  Richard  Vogel, graduate  student  at  Caltech. 

Mr.  Carl  Kenreich,  graduate  student  at  Caltech. 

Mr.  K.  Young,  student  at  U.C.L.A. 


RAYLEIGH  WAVE  TRANSMISSION  ON  ELASTIC  WEDGES* 


J.KAXEt  AND  J  .  S  I*  E  X  C  K  f 


A  method  is  introduced  liv  nhivh  transmission  nr  Conversion  coefficient*  can  readily  l>c  calculated  for  surface-wave 
diffraction  at  crustal  >li.  Continuities  which  can  be  idealized  as  wcdgi.-shaj>cd  obstacles.  The  method,  an  iteration 
procedure,  is  illustrated  by  an  application  to  Rayleigh-wave  propagation  on  elastic  wedges  with  free  lvoundaric*. 
It  is  sho" ii  that  ihe  leading  term  of  the  iteration  j*r  --cdure  yields  an  espressn  :  for  the  transmission  coefficient 
w  hich  matches  lhec\|teri,iieiilal  i.ala  for  a  range  of  \u  .ge  angles  w  hich  includes  the  crustal  wedges  commonly  found 
in  seismic  practice,  Numencal  data  are  presented  whh  n  illustrate  the  transmission  coefficient's  variation  in  terms  of 
Poisson  s  ra  i  io  and  the  w  t  dgc  angle. 


INTRODUCTION 

The  diffraction  of  Rayleigh  waves  by  w  riges  is 
a  topic  of  practical  importance  because  it  is  re¬ 
lated  to  the  generation  of  seismic  noise  and  to 
other  geophysical  and  acoustical  applications. 
In  addition,  it  has  considerable  inherent  theoreti¬ 
cal  interest  because  the  associated  boundary 
value  problem  has  not  yet  yielded  to  any  avail¬ 
able  mathematical  techniques.  In  this  report,  we 
introduce  a  procedure  which  is  useful  for  estimat¬ 
ing  the  propagation  effects  on  seismic  waves  prop¬ 
agating  past  wedge-shaped  obstacles.  Our  method 
is  basically  an  iteration  procedure  by  which  the 
boundary  conditions  on  either  face  of  the  wedge 
are  alternately  satisfied.  The  convergence  of  the 
associated  series  and  the  higlicr-order  iterations 
will  be  discussed  in  a  subsequent  paper.  It  is  the 
purpose  of  this  report  to  demonstrate  that  .he  first 
iteration  yields  the  principal  information  desired 
by  the  seismologist,  namely  the  field  that  propa¬ 
gates  past  thediscor.iir.uity.  This  data  is  obtained 
readily  by  our  a:valj  .use  live  faces  of  wedges 
associated  with  crus;.,,  i.. -continuities  arc  nearly 
parallel,  and  our  procedure  yields  the  transmis¬ 
sion  c  .  civnts  by  rather  elementary  operations. 
Wc  shah  illustrate  the  treatment  by  calculating 
the  transmission  coefficient  of  a  Rayleigh  wave 
which  is  incident  upon  a  free  elastic  wedge. 

In  sonic  respects  our  first  iteration  is  related  to 


the  Kirchoff-Huygcn  method  which  finds  an 
approximate  field  by  integrating  a  suitable 
Green's  function  over  an  aperture  with  a  weight 
function  proportional  to  an  assumed  illumination. 
Our  method  differs  iu  several  respects:  wc  neither 
need  obtain  a  Green's  function,  nor  integrate  it, 
and  our  calculations  can  be  systematically  im¬ 
proved. 

roaxuLATioN 

We  consider  an  isolated  wedge  with  free 
boundaries  (Figure  1).  It  is  very  convenient  to 
introduce  two  coordinate  systems:  first,  the 
x,  y  system  simplifies  the  discussion  of  the  left 
face,  and  second,  the  {,  ij  system  is  appropriate 
for  the  right  one.  The  two  coordinate  systems  dif¬ 
fer  by  a  rotation  through  an  angle  6  which  we  shall 
call  the  wedge  angle',  thus,  wc  have 

x  =  {  cos  9  —  ij  sin  8, 

(1) 

y  -  £  sin  8  +  ij  cos  8, 

and  the  inverse  transformation 

£  "  x  cos  9  +  y  sin  8, 

(2) 

ij  “  —  x  sin  0  +  y  cos  9. 

The  angle  8  is  positive  when  measured  counter¬ 
clockwise  from  $**  0,  the  x  axis,  so  that  negative 
values  cf  8  correspond  to  clastic  wedges  whose 


*  Thc  .rc.5.l'jrch  described  in  this  report  ha*  been  sponsored  by  the  Advanced  Research  Projects  Agency,  Project 
Vela-Uniform,  ur.ner  contract  AF  19(628)419  o i  the  Gcophy-ics  Research  Directorate,  Air  Force  Cambridge 
Research  Lab u--:,- ,:iei.  Office  of  Aerospace  Research,  United  State*  Air  Force,  Bedford,  Massachusetts.  Matsu- 
script  received  rive  Editor  January  25, 1963. 

f  University  of  ....ode  Island,  Kingston,  R.  I. 
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_  ,  ,  „rol>!cm-  A  Rayleigh  wave  L  incident  from  *-  -  -  along  the  left  face  of  a  free  elaatre 

Fic.  1.  The  geometry  of  jn  thisdhistralion,  the  wedge  angle  #  is  negative. 


interior  angles  are  less  thane  radians  and  positive 
values  to  wedges  whose  interior  angles  arc  greater 
than  e  radians.  If  «- 0,  then  the  wedge  degen¬ 
erates  into  a  half-space  with  no  discontinuity  at 
the  origin. 

In  the  absence  of  body  forces,  small  displace¬ 
ments  l(u,  r,  w)  of  an  elastic  solid  characterized 
by  the  Lame  parameters  X,  m.  and  density  p,  can 
be  derived  fro  a  a  scalar  potential  t>,  and  a  vector 
potential  V(^„^S,W,  namely, 

g(«,  r,  w)  «  grad  *  +  curW’^i,  h,  *»)• 

For  tw.  dimensional  motions,  $«$(*,  y), 

-  ;  v  ),  and  we  can  neglect  pure  distortions 
.  v  littinE  tf-i-tfs-O.  *o  that  without  introducing 
ny  confusion  we  can  drop  the  subscript  on 
j,.(.r  v).  For  monochromatic  vibrations,  we  can 
suppress  a  time  factor  t~~',  and  it  can  be  shown 
that  *  and  *  satisfy  the  reduced  wave  equations; 


(V5  +  V)*  -  0,  *«5  -  «V'(X  +  m),(3) 

(V*  4-  !■;■)■;  =  ot 

i.i  the  ir-'e.-.or  ot  the  wedge.  Let  '0  be  dummy 
variables  which  represent  either  the  (£,  >i)  or 
(x?  y)  coordinate  systems;  then  the  tangential 

«r«*i  ... 


and  the  normal  stress, 

(d'4>  av  \ 

=  +  J.  O 

must  vanish  on  either  face  of  the  wedge.  As  our 
excitation,  wc  shall  choose  a  Rayleigh  wave  of 
unit  amplitude, 

R(4>r,  i.)  *  +  +r, 

in  the  elastic  wedge  which  is  incident  from 
x  —  —  " ,  so  that 

i*t  = 

*(♦„  ir)  -  \  (  ) 

=  2c'"4''*’  *•  *• 

Wc  shall  speak  of  the  coefficient,  in  this  case 
unity,  of  the  comprcssional  potential  <?r  as  the 
amplitude  of  the  Rayleigh  wave  The  parameters 
kr  and  Z  are  so  chosen  that  the  Rayleigh  wave 
exerts  neither  normal  nor  tangential  stress  at 
the  free  clastic  surface  at  v**0.  In  general, 
j  JLj  >  |  k.\  >  |  k,\ ,  so  that  both  «r  and  ^  are  ex¬ 
ponentially  attenuated  with  increasing  depth. 
We  seek  diffracted  fields  (&,  W  which  wc  ca" 
add  to  the  Rayleigh  wave  so  that  the  total  held 

ur,  where 


“  ’  (  aV  d'+*H\  “r,  where 

p»{4>,  t)  -  f*  y  dldl,  ~  a,»  +  a„V  ’  «r  -  (*-  +  *<)  +  W*  + 
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exerts  neither  normal  nor  tangential  stress  on 
either  wedge  face. 

THR  ITtR/TCD  SOI.UTION 
Preliminary  remarks 

Let  «|W'  denote  the  incident  held,  or  zeroth- 
ordcr  iteration,  which  in  our  ease  would  be  the 
Rayleigh  wave 

«!»«'  ®*  <t>r  +  'i'r- 

Each  term  in  this  expression  is  a  wave  function, 
and  the  combination  4>r-\-\pr  satisfies  the  bound¬ 
ary  conditions  on  the  left  face  of  the  wedge,  but 
not  those  on  the  right  face.  There  is  no  known1 
way  by  which  one  can  add  wave  functions 
to  *!„•  so  that  all  boundary  conditions  are  satis¬ 
fied.  However,  it  is  possible  to  introduce  wave 
functions  (fr-.tf'1)  so  that 

-  «lo.*  +  (*',  *')  =  (*,  +  -f  *•) 

satisfies  the  boundary  conditions  on  the  right  face. 
One  can  think  of  0'  and  <fil  as  elastic  waves  gen¬ 
erated  by  the  response  of  the  second  face  to  the 
incident  Rayleigh  wave.  Of  cou.se  u'  will  not 
satisfy  the  boundary  condition  on  the  left  face 
because  the  terms  (*',  V)  have  been  introduced 
without  consideration  of  th;  left-hand  boundary 
conditions.  However,  on  heuristic  principles,  we 
would  imagine  that  if  the  wedge  angle  8  were 
small,  then  very  little  of  the  secondary  energy 
associated  with  (*',  *')  would  be  returned  to  the 
left  lace.  Thus,  although  we  could  introduce  an¬ 
other  term  in  the  iteration  procedure,  namely 

si’.*=  k*  -f  (if,1,  f1), 

so  that  «*  satisfies  the  left  boundary  condition, 
we  shall  not  do  so  in  the  present  paper,  beca  ise  we 
shall  show  that «'  yields  a  transmission  coefficient 
which  agrees  very  well  with  experimental  data 
for  a  rather  wide  variation  of  8. 


Pls(4>„  f,)  /($),  (7) 

,  ir)  =*  g({),  (8) 

to  that  half-plane  ij-O,  {£0.  In  our  context,  the 
expressions  p (,  and  are  linear  differential 
operators  on  if,  and  ^  which  involve  exponential 
functions  of  algcbraicargunent.  Since  differentia¬ 
tion  of  such  exponentials  is  equivalent  to  multi¬ 
plication  by  algebraic  expressions,  it  is  very  con¬ 
venient  to  introduce  the  following  functions  of  k 
defined  over  the  complex  wave  number  plane: 


P(k) 

P. 

(9) 

Q(k) 

-  Mk'  -  *.«), 

(10) 

m 

-  -  X*.*  -  -  *«), 

(ID 

S(k) 

«  -  2knVkJ~- 

(12) 

and  a  useful 

determinant 

Del 

,«)- 

Q(k)  I 

(13) 

*(*) 

s(k)  r 

It  is  important  to  note  that  since  the  Rayleigh 
wave  /?($»,  \pr)  exerts  neither  normal  nor  tan¬ 
gential  stress  at  y-0,  we  have 


Pik,)  +  Q(kr)Z  -  0,  (|4) 

R(k,)  -f  *{kr)Z  -  0,  (15) 

so  that  the  columns  of  Det(*,)  are  proportional, 
or  the  period  equation 

Det  (*)  *  0  (16) 

has  a  root  at  it-*,. 

With  these  definitions,  the  Rayleigh  wave 
stresses  /(()  and  g(()  on  the  right  face.  g— 0, 
f  “0,  are  explicitly 

/({)  ~  P(a)e««  +  Q(0)Ze*t  (17) 
and 


Integral  representations 

The  incident  Rayleigh  wave  u,M  does  not  sat¬ 
isfy  the  boundary  conditions  on  the  second  wedge 
face  for  it  applies  nonvanishing  stresses: 

*  Indeed,  because  of  the  geometric  singularity  at  the 
sharp  edije,  there  may  well  be  no  solution  if  the  bound¬ 
ary  conditions  are  used  uniformly  at  all  distances  from 
the  vertex.  The  only  well-posed  boundary  value  problem 
of  this  type  may  be  one  where  the  w  edge  has  a  rounded 
tip- 


l({)  -  R(a)e“l  -f  S(0)Ze**,  (t8> 

where  a(8)  and  0(8)  lie  on  ellipses  (see  Figure  2) 
in  the  complex  k  plane: 

a(8)  *=  kr  cos  8  -  Wk?  ~  *.*  sin  8,  (19) 

fi(8)  =  kr  cos  8  -  tvV  -  *,*  sin  8.  (20) 

Xote  that  for  small  values  of  8,  both  «  and  0  are 
approximately  equal  to  kr.  For  {>n,  we  can  re¬ 
write  the  streoses  (17)  and  (18)  in  integral  form 


*ti  i  wi  fwrtffrrrtnrcrwi™  mtifflpnfn,mimnnrt 
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by  the  Cauchy  integral  theorem  as 


under  the  integral  sign 


I  C  (  A»  S{ti)  \ 

where  the  contour  C  is  sketched  in  Figure  3.  Ex¬ 
pressions  (21)  and  (22)  display  the  stresses  on 
q»0  as  Fourier  transforms.  Let  us  note  that  if 
the  boundary  value*  «({,  0)  of  a  function  u((,  q) 
are  given  i>i  the  form 

u((,  0)  ■  f  U(k)e<*Jk,  I  >  0, 
2wiJ  c 

then  it  is  a  trivial  matter  to  analytically  continue 
tt((,  0)  at  a  solution  of  the  wave  equation 

(V*  +  «')«({,  n)  -  0 

into  the  half-space  q  <0  by  inserting  a  factor 


and  suitably  indenting  ’he  contour  C  about  the 
branch  points  ±*  in  the  complex  k  plane.  Vertical 
branch  cuts  have  been  chosen  so  that  the  cui  k 
plane  (see  Figure  3)  represent*  that  sheet  of  the 
Ricnunn  surface  for  which  the  radicals  \'k,,-kt 
and  i'are  respectively  +i,and  -ft,  at  the 

origin  F»0;  by  analytic  continuation,  this  fixes 
the  values  of  the  radicals  throughout  the  entire 
cut  k  plane. 

With  these  preliminaries,  we  seek  the  first  un¬ 
known  wave  functions  ($\  /-*)  in  the  form: 

-  --  f  (23) 

2  nJ  c 

m^-^Mk,  (24) 

2 ttJ  c 

where  the  kernels  X(k)  and  Y(k)  are  to  be  chosen 


j3(0) 


k?-k2, 


-> 

Re  k 


Fic  2  The  location  of  the  poles  *-a.  3  in  the  complex  *  plane  as  a  function  of  the  wedge  angle*.  Becaw  of  the 
minus  signs  in  equations  (19)  and  (20),  the  poles  a  and  fl  he  in  the  upper  half  plane  when  i  »  negause. 
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to  that  the  combination  1>f+V)  exerts 

neither  Ungen  rial  stress 

fc,($r  +  4>\  i'  +  V)  “  0,  (25) 

nor  normal  stress 

r  +  ir  +  i')  -  0  (26) 

on  tin  wedge  face  ij-0,  {>0.  Assuming  that  we 
can  uitlerentiate  freely  under  the  integral  signs, 
we  find  that  (25)  and  (26)  require  that  -Y(i)  and 
Y(k)  be  solutions  of  the  pair  of  algebraic  equa¬ 
tions: 


P(k)X(k)  +  Q{k)Y{k) 

m  r  pm 

“  ”  L*  - « 

/?(*)*(*) +J(*)K(4) 

fX(a) 

"  ~  U  -  « 


+  Z 


+  z 


QW 1 

k-er 


1 

k-el' 


(27) 


which  can  be  wived  by  inspection  as 


*(M ) 

Y(k,e) 


A(k,9)S{k)  -  B(k,6)Q(k) 
Del  (*) 

B(k,e)p(k)  -  ■-!(*,  e)R(k) 
Det  (4) 


»  (28) 
*  (29) 


where  the  9  variation  is  conUined  in  the  expres¬ 
sions 


t.id 


A{k,  9) 

B(k,  6) 


PM 

k  —  a 

m_ 

k  —  a 


+  Z 


Z 


Q(&) 

k-fi' 

m 

k-0 


(30) 

(31) 


Hence,  the  field  up  to  the  first  iteration  would  be 


x  (Plant  tun  the  pH*  ft'  tqutlitH  32.) 
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«•(♦'.  +')  “ 


*■ + ivJ, 

>-+u 


A>.W>  - 


2 wiJ  c  (*) 

1  f  B(k,8)P(k)  -  A{k,$)R{k) 
2 wiJc  Det  ( k ) 
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C9i*y<irison  with  experiment 

The  integrals  in  (32)  represent  the  first-order 
reaction  of  the  right  wedge  face  tn  the  incident 
Rayleigh  wave,  and  include  both  body  waves  and 
surface  waves  as  contributions  to  the  net  field. 
The  cylindrical  body  waves  can  be  evaluated  by 
a  standard  appeal  to  the  method  of  saddle  points 
and  will  not  be  discussed.  The  surface  waves 
emerge  as  residue  fields  of  the  various  poles  of 
(32).  By  construction,  th:  poles  at  i-a  and 
give  rise  to  a  field  which  just  cancels  the 
incident  Rayleigh  wave  on  the  right  face  of  the 
wedge.  We  have  observed  in  equation  (16)  that 
DET(t)  has  a  zero  at  k  -  i„  so  that  both  integrals 
in  (32)  have  another  pole  a;  4  “  if.  Clearly ,  the 
residue  field  associated  with  this  pole: 


-  T(6,  <r)<rK,<v'*’'-‘«\ 

_ (33) 

=  T(B, 

represents  the  leading  contribution  to  the  trans¬ 
mitted  Rayleigh  wave.  The  amplitude  coefficient 
T{0,  a)  of  this  transmitted  field  can  be  readily 
calculated  from  (32)  as 

A{k„  0)5(4,)  -  B(kr,eMkr) 

T(S,  <r)  - - - j - ~  ’  (34> 

_  Det  (ft) 
dk  Is-*, 

which  is  i  function  of  the  wedge  angle  0  and 
Poisson’s  ratio  <r.  When  there  is  no  discontinuity, 
it  can  easily  be  shown  that 

lim  T(9,  a)  *  1,  (35) 

•-« 


.  . .  Knopoff  end  Gang! 
—  Calculated  Values 


•120’  -ICC  *80“  -W  -4 <f  -20’ 


*»  — • 


100*  120’  140'  160*  ISO* 


fie.  4.  a  .i  k*c  m>  i".  ■"  “• 

of  Knopoff  and  Gftnp  for  Poisson  i  ratio 
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Fio.  3.  The  absolute  square  of  the  transmission  coefficient  T,  which  a  an  even  function  of  the  wedge  angle  t. 

^as  it  should.  Furthermore,  since  both  the  num-  0  >  +  70®,  some  of  the  discrepancy  should  be  re- 

erator  and  denominator  of  (34)  are  homogeneous  moved  by  including  the  contribution  of  higher 

functions  of  the  same  order  of  the  parameter  kr,  iterations,  but  experimental  sources  of  error  which 

there  is  no  frequency  dependence  in  the  leading  increase  the  transmission  coefficient  should  not  be 

term  of  the  transmission  coefficient.  ruled  out:  these  might  include  a  slight  rounding 

The  range  of  utility  of  (34)  can  be  gauged  by  of  the  vertex,  or  multiple  scattering  from  minor 

an  appeal  to  carefully  measured  experimental  inhomogeneities  in  the  model.  A  much  more 

evidence:  in  Figure  4,  we  plot  |  T| 1  as  a  function  interesting  item  in  the  experimental  data  is  the 

oiO  for  Poisson's  ratio  a  -0.266,  and  compare  it  to  secondary  peak  in  the  transmission  coefficient  at 

the  data  of  Knopoff  and  Gangi  (I960).  We  note  about  0»  — 100*  which  is  not  reproduced  by  (34), 

that  the  experimental  points  straddle  the  calcu-  the  first-order  transmission  coefficient.  Similar 

lated  values  for  —60*  <6  <+70®.  Oi  course,  this  humps  appear  even  in  less  refined  measurements 

range  includes  practically  all  crustal  wedges  com-  (deBremaecker,  1958;  Viktorov,  1958),  which 

monly  found  in  seismic  practice.  implies  that  this  secondary  maximum  is  a  charac- 

However,  there  are  theoretical  reasons  why  ex-  teristic  feature  of  Rayleigh-wave  diffraction  by 

treme  values  of#  are  of  interest,  and  it  maybeuse-  wedges.  Although  this  peak  docs  not  appear  in 

ful  io  comment  upon  Figure  4  in  more  detail.  For  our  first-order  analysis,  a  qualitative  inspection  of 


. . . 
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F;o.  6  The  phase  of  T(t,  <r)  which  is  an  edd  function  of  the  wedge  angle  9.  Positive  values  of 
phase  correspond  to  positive  values  of  I. 


the  higher-order  iterations  indicates  a  theoretical 
basis  for  this  phenomenon  in  terms  of  the  sensitive 
dependence  of  the  transmission  coefficient  upon 
the  location  of  the  reflected  Rayleigh-wave  pole 
when  —  11O*<0<  —  80*.  This  invites  an  inter¬ 
pretation  in  terms  of  a  resonance  between  the 
reflected  and  transmitter!  fields  within  this  range 
of  9. 


Prescnlolion  of  numerical  data 

The  agreement  between  theory  and  experiment 
for  the  range  —  60*  <0<-t-7O°  prompts  a  de¬ 
tailed  display  of  the  behavior  of  the  transmis¬ 
sion  coefficient  T(9,  o)  in  this  band.  In  Figure  5, 
we  illustrate  the  dependence  of  |  T(8,  ff)| which 
is  an  even  function  of  6,  upon  Poisson’s  ratio  <r. 
For  fixed  9,  the  magnitude  of  the  transmission 
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coefficient  varies  inversely  a»  Poisson’s  ratio. 
However,  we  note  that  this  variation  is  not  pro¬ 
nounced.  In  general,  the  transmission  coefficient 
is  a  complex  number,  and  in  Figure  6  we  plot  the 
phase  of  T{8,  a),  which  is  an  odd  function  of  8 
for  the  same  range  of  0  as  in  Figure  5.  This  phase 
variation  is  altnqst  linear  and  is  practically  in- 
sensitive  to  changes  in  Poisson's  ratio. 

CONCLUSIONS 

The  method  vee  have  used  to  discuss  Fiylcigh- 
wavc  transmission  can  be  readily  adapted  to  esti¬ 
mate  transmission  or  conversion  coefficients  for 
surface-wave  diffraction  at  other  crustal  discon¬ 
tinuities  which  can  be  :,,cali*ed  as  wedge-shaped 
obstacles.  The  principal  advanttge  of  the  tech¬ 
nique  is  that  these  coefficients  can  be  obtained 


with  an  economy  of  effort,  particularly  if  the 
wedge  angle  is  small.  Although  contour  integra¬ 
tion  has  been  Introduced  In  describing  the 
method,  an  examination  of  the  analysis  will 
reveal  that  the  residue  coefficients  could  have 
been  calculated  directly  from  the  boundary  con¬ 
ditions,  and  by  an  explicit  representation  of  the 
incident  field. 
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The  Theory  of  Surface  Wave  Diffraction  by 
Symmetric  Crustal  Discontinuities 
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Summary 

The  techniques  of  functional  analysis  arc  used  to  find  approximate 
solutions  to  the  p'oblem  of  surface  wave  diffraction  by  symmetric 
wedge-shaped  obstacles.  The  complex  amplitudes  of  the  scattered  surface 
waves  arc  found  by  a  variational  procedure  which  requires  but  the 
solution  of  a  small  number  of  linear  algebraic  equations:  the  Gauss 
normal  equations.  The  analysis  is  illustrated  by  a  discussion  o  ' the  vector 
problem  of  Rayleigh  wave  diffraction  by  clastic  wedges,  and  the  scalar 
problem  of  Love  wave  diffraction  in  a  layered  wedge  which  supports 
two  propagaL  ics. 

Introduction 

A  major  barrier  in  comparing  seismic  t'neo  ,  with  observed  wave  trains,  stems 
from  the  fact  that  elastic  wave  characteristics  arc  influenced  by  discontinuities 
along  the  propagation  path,  and  any  understanding  of  such  sign.-'  corruption  would 
require  a  knowledge  of  the  diffracted  fields  at  the  appropriate  .aclcs.  However, 
even  a  relatively  sir  . pic  crustal  feature  such  as  a  discontinuous  change  in  terrain 
presents  major  difficulties  if  the  relevant  problem  of  diffraction  in  a  wedge-shaped 
region  is  considered  Although  some  first-order  calculations  have  been  made 
by  Lapwood  (1961),  Kane  &  Spence  (1963),  and  Hudson  <!i.  Knopoff  (1963),  the 
theoretical  discussion  of  the  diffraction  effects  are  hampered  by  the  intractability 
of  the  associated  boundary  value  problems.  In  this  report,  we  'l.ovv  how  one 
can  take  advantage  of  symmetry  considerations  and  variational  techniques  to 
estimate  rapidly  reflection,  transmission,  and  conversion  coefficients  for  clastic 
wave  diffraction  at  symmetric  vvedec-shaped  obstacles.  In  Part  1,  we  illustrate  the 
ideas  by  a  discussion  of  the  vector  pioblcm  of  Rayleigh  wave  propagation  along 
the  faces  of  an  elastic  wedge  "  ith  free  boundaries.  In  Part  2,  we  analyse  the  scalar 
problem  of  multi-mode  Love  wave  diffraction  in  a  symmetric  layered  wedge. 

Part  1.  Rayleigh  waves  on  an  elastic  wedge 
1.  fundamental  equations 

The  tr  mors  s(ii,  v,  w)  of  an  clastic  solid  characterized  by  the  Lame  parameters 
/,  pi,  and  density  p,  can  be  derived  from  a  sealer  potential  <}>(x,y,z),  and  a  vector 
potential  V(i/\(.v.  y,  ;),  >’■  vMv>  L  ')J  by  the  relation 

s(«,i,w)  =  V.tf+Vx'F.  (LI) 

*  Received  in  oritin*!  form  t°63  November  8. 
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For  iwo-dimensional  motions  which  arc  independent  of  the  : -coordinate,  both 
<t>  and  S’  are  but  functions  of  v  and  y,  or  f/>  =-  4>{,\,  j  ).  and  S'  ~  S'(.t,  >•),  Furthermore, 
we  can  neglect  pure  distortions  by  setting  ipt  <=  *=  0,  so  that  the  vector  potential 

S'  -  S’[0. 0.  ip,(x,  >•))  is  characterised  by  one  scalar  component  and  the  subscript 
on  ifr,  can  dropped  without  confusion.  If  we  assume  that  the  vibrations  are 
harmonic,  we  can  suppress  a  time  factor  and  it  can  be  shown  that  ^  and 
the  ^-component  of  the  vector  potential,  satisfy  the  reduced  wave  equations 

+  =  0,  k}  =  vi'pftk+nl  (1.2) 

(V,+fc,I)v!'  -  0.  k}  =  u)lpln .  (1,3) 

Once  the  potentials  </>  and  1 1>  arc  known,  the  displacement  vector  s  is  given  by  (1.1), 
and  the  normal  and  tangential  stresses  a„,  o,}  can  be  obtained  from  the  stress 
dyadic  £($,  ^)  given  in  symbolic  notation  as 

<c{<t>,4>)  =  ^3V.§  +  p(Vs  +  sV)  (1.4) 

where  3  is  the  unity  dyadic,  the  idcmfactor. 

2.  Hie  Rayleigh  wave  potentials 

A  time  harmonic  Rayleigh  wave,  or  if- wave  for  brevity,  is  comprised  of  a  pair  of 
exponential  solutions  of  (1.2)  and  (1.3),  j?  =  3f(0„,  tfr*).  If  the  elastic  solid  lies 
within  the  half-space  y  £  C  ^  ;  Figure  1(a)),  then  these  solutions,  in  polar  coordin¬ 
ates  x  »  rcosS,  y  «  rsinfl,  assume  the  form 

'$*  -  «***  (1.5) 

-  -iT****  (1.6) 

wherein  the  sponential  variation  is  given  by 

0(3)  =  -(v/JV')1 1  sin  3  +  i cos  3,  v}  «•  (1.7) 

P 

AS)  =  VH  “(fa/v.)*]  sin  3  -»• » cos  3.  vf  *  -  (1.3) 

P 

and  T,  the  magnitude  of  the  ratio  of  the  shear  potential  coefficient  of  to  the 
comprcssional  one  is 

i ( ,9) 

and  E  is  a  dimensionless  dictance  parameter 

•=•  *=  kKr  =  -r.  (1.10) 

o* 

The  parameters  t’{  and  v,  represent  the  velocities  of  the  conp.essional  and  shear 
body  waves  respectively.  For  a  given  Poisson’s  ratio  a 


2{X+p) 

one  needs  solve  for  the  Rayleigh  wave  velocity  t„,  which  is  less  than  rt.  so  that  the 
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stresses  induced  by  </>*  and  yR  vanish  at  the  surface  y  -  0 


o,M*.  0*>  =  o,A<t>ie  0*)  =  0,  y  =■  0. 


(M2) 


With  this  choice  of  i*.  the  </>„  and  given  by  (1.5)  and  (l/»)  arc  the  vector  and 
scalar  potentials  characterizing  a  Rayleigh  wave  travelling  to  the  right  with  unit 
amplitude.  (We  shall  stvak  of  the  coeDicient  of  the  compressional  potential  as  the 
amplitude.)  It  is  very  convenient  to  note  that  we  can  reverse  the  direction  of  uny 
hartnonic  wave  by  the  operation  of  complex  conjugation.  Thus 


(  R4>'n  -  Re''™ 

t/'*)  =  I 

{  R&'k  -fiRlVv,,s 


(1.13) 

(114) 


represents  a  Rayleigh  wave  travelling  in  the  opposite  direction  with  amplitude  R. 


2  Formulation  of  the  boundary  value  problems 


(a)  The  major  problem 

The  conundrum  posed  by  Rayleigh  wave  diffraction  in  a  wedge-shaped  region 
is  to  find  additional  solutions  and  ^  of  (1.2)  and  (1.3)  >  hich  represent  diffracted 
fields  i*i  the  interior  of  the  wedge  such  that  the  stresses  <r„  and  <rr3  vanish  on  both 
faces  of  the  wedge  for  Jt- wave  excitation  along  one  face.  In  our  geometry,  (see 
Figure  1(a)),  the  df-wavc  is  incident  from  infinity  along  the  negative  x-axis.  We 
shall  be  principally  concerned  with  calculating  the  complex  amplitudes  of  the  re¬ 
flected  and  transmitted  JP-waves  as  a  function  of  the  wedge  angle  0  and  Poisson’s 
ratio  a.  This  task  can  be  substantially  eased  by  reducing  the  major  problem  to  a 
pair  of  minor  problems  involving  even  and  odd  symmetric*  Consider  the  incident 
£f-wavc  along  the  left  wedge  face  to  be  the  sum  of  two  waves,  each  of  half  amplitude. 
Likewise,  the  absence  of  any  excitation  along  the  right  wedge  face  is  equivalent  to 
a  pair  of  incident  Jf- waves  along  it,  each  of  half  amplitude,  but  opposed  in  sign. 


Fig.  I. — A  unil  Rayleigh  Wave  incident  along  one  face  of  an  elastic  wedge 
is  equivalent  to  four  partial  waves.  The  partial  waves  of  like  parity  comprise 
the  excitation  for  the  even  problem.  a-sJ  the  partial  waves  of  unlike  parity 
furnish  the  initial  disturbance  for  the  odd  problem. 
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(see,  Figure  1(h)).  The  four  partial  ^ -waves  can  be  separated  into  two  groups. 
First,  a  pair  of  Jf-waves  on  either  face  of  like  parity  which  serves  as  the  excitation 
of  what  we  call  the  even  minor  problem.  Second,  another  pair  of  Jf-waves  whose 
amplitudes  are  of  unlike  parity  which  constitutes  the  excitation  of  the  odd  minor 
probknL  If  we  designate  the  diffracted  fields  of  the  even  and  odd  problems  by  the 
subscripts  e  and  o  respectively,  then  since  we  arc  dealing  with  linear  equations,  the 
desired  major  potentials  <f>4  and  can  be  expressed  as  a  superposition  of  the 
minor  potentials 

=  +  (1. 15) 

-  tf'.-Mf'. .  .  (1.16) 

and  likewise  the  displacement  vector  s  —  s,  +s„  can  be  decomposed  into  even  and 
odd  components. 

(6)  The  even  minor  problem 

In  the  even  problem,  the  wedge  will  suffer  only  ever,  displacements*  s,  about  the 
plane  of  symmetry,  and  as  a  result,  there  can  be  no  component  of  normal  displace* 
ment  along  the  plane  of  symmetry  at  3  =  0/2 -a.  It  follows  that  the  even  problem 
is  equivalent  to  finding  the  potentials  in  a  bisected  wedge  with  one  face  free  of 
stresses  which  supports  the  incident  Rayleigh  wave,  and  the  other  face  so  con¬ 
strained  that  the  normal  displacement  vanishes  there.  That  is,  we  seek  solutions 
&  and  of  (1.2)  and  (1.3'  ;n  a  wedge  of  half-angle  0/2 

'<U0*  +  0., -  ^/H-0„<i>  +  ^,)  -  0,3  -  it  (1.17) 
even  j  d  ^ 

+  +  |;W*  -M^.)  «  0,  3  -  0/2-  a 


(c)  The  odd  minor  problem 


(118) 


By  similar  arguments,  the  odd  problem  which  involves  $r  is  a  complemen¬ 
tary  version  in  the  bisected  wedge,  wherein  the  tangentir  I  displacements  must  vanish 
identically  along  the  plane  of  symmetry,  i.c., 


+  +  »=  0,  3-11  (1.19) 

d  I  c 

^W*  +  0O)+  +  =  0,  3  »  0/2— a  (1.20) 

(d)  The  reflection  coefficients 

For  either  the  even  or  the  odd  problem,  the  solution  will  contain  a  reflected 
Rayleigh  wave.  Let  p,(0/2,  o)  and  pQ(0/2,  o)  be  the  complex  rat  •$  of  the  reflected 
to  the  incident  Rayleigh  wave  amplitude  for  the  even  and  odd  minor  problems  in  the 
bisected  wedge.  The  reflection  coefficient  R(Q, o)  for  the  original  major  problem 
will  be 

R(Q.o)  =  i(p,(0/2,  a)  +  po{0/2,  a)],  (1.21) 

and  likewise  the  overall  transmission  coefficient  T(0,  o)  will  be 

ne.  O)  =  i(p,(0/2, ,-)  -  p,(0/ 2,  er) J.  (1.22) 

Tlie  displacements  and  lhe  compressions!  pnicnual  will  be  even  about  iiie  plane  of  jymmeiry,  but 
tht  ihwf  potential  well  be  an  odd  function,  and  rice  versa  for  the  odd  problem. 


th<  Ihrun  of  writer  »sit  ditTraeiton 


Formulas  (1.21)  and  (1.22)  can  be  verified  by  a  glance  at  gurc  ((/>)  which  indicates 
that  the  overall  reflection  coeflicicnt  R  result'  from  a  superposition  of  the  partial 
reflection  coefficients  }(/>,  +  />„),  and  the  transmission  coefficient  Tfrom  their  inter¬ 
ference  }(/>,-  pv). 

4.  The  variational  principle 
(a)  Discussion 

Variational  procedures  consist  of  assuming  a  suitable  trial  function  containing 
unspecified  coefficients,  and  then  choosing  these  parameters  to  minimise  certain 
quantities.  One  major  advantage  of  the  variational  method  is  that  first-order 
accuracy  in  the  trial  function  usually  gives  results  which  arc  accurate  to  second- 
order,  because  of  the  stationary  character  of  the  approximation. 

A  natural  aperture  in  the  present  problem  is  the  plane  of  symmetry  and  we  can 
assume  it  to  be  illuminated  by  an  incident  Jf-wave,  and  a  reflected  one  with  an 
adjustable  amplitude.  In  the  even  problem,  the  net  angular  displacement  must 
vanish  along  the  plane  of  symmetr,.  A  unit  J’-wave  travelling  to  the  right  gives  rise 
to  the  angular  component  of  the  displacement 


and  likewise  an  ^-wave  of  amplitude  p,  moving  to  the  left  generates  the  disturbance 

PJ  t*  =  (1.24) 

which  apart  from  an  amplitude  factor  is  the  complex  conjugate  of  (1.23).  Oni  >  if 
there  is  no  discontinuity,  or  if  0  =  a  can  we  make  the  angular  displacement  of  the 
trial  function 

il  =  slT+P.WT  •  (1.25) 

vanish  for  all  r  along  the  plane  of  symmetry  by  properly  choosing  pr  Otherwise 
a  and  /?  are  complex  fl*.  and  no  choice  of  pt  car.  make  s\  vanish  at  more 

•ban  an  isolated  set  of  points.  There  arc  at  least  two  ways  by  which  we  can  improve 
matters.  We  could  use  a  more  complex  trial  function  which  acknowledges  body- 
wave  contributions  to  the  diffracted  field,  or,  since  the  residua!  displacement  sj  is 
explicitly  known,  we  can  use  it  as  the  aperture  illumination  or  a  Green’s  theorem 
type  calculation  to  correct  the  variational  estimate. 

Howtver,  the  practical  seismic  interest  is  in  the  realm  of  small  discontinuities, 
and  for  this  case  we  shall  see  that  the  elementary  trial  function  yields  satisfactory 
results. 

(b)  Definition  of  the  scalar  product 

While  there  are  many  gauges  by  which  sj  can  be  minimized,  we  shall  choose  a 
p,  which  minimizes  s]  in  the  mean  square  sense.  For  this  purpose,  let  us  define 
the  complex  scalar  product  of  two  functions  u(r,3)  and  v(r,  3)  to  be 

r 

(u,  y)  =  J  u(r,o/)i,(r,3)dr,  (1,26) 

where  the  integration  is  to  be  carried  out  along  the  plane  of  symmetry.  3  =  0/2 -n. 
To  each  complex  function  u(r,9)  we  can  attach  a  positive  definite  number,  the 
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norm  of  u  or  jjujj  which  is  defined  to  he  v  («.  a).  The  norm  fiijjd  pends  on  th  wedge 
angle,  and  is  to  be  distinguished  from  u1  -  (u.ir)  which  is  in  general  complex. 


(c)  The  even  subsidiary  reflection  coefficient 

With  this  notation,  the  mean  square  value  of  the  angular  displacement  of  the 
trial  function  sj  is 


KT  -  tf'+pjtirr.fr+pjiJm 

and  this  will  be  a  minimum  if,  and  only  if,  p,  is  chosen  as 

,  (j'n‘  . 

rjei2.nl -jjpp. 

or  explicitly  in  terms  of  a{9),  /?(S)  and  T, 


P,(Q/l°) 


1  /da'' 

\*  2iVpj 

“2«V*9y 

I  +  2  a+  \ 

_LJ 

fcTrW 

+irf 

V  (M\ 

«+ a*j 

d9  P+F 

*+F\d») 

(1.27) 

(1.28) 


(129) 


(</)  The  odd  subsidiary  reflection  coefficient 

Jf  we  use  an  analogous  trial  function,  and  similar  reasoning,  we  find  that  if  p .  is  to 
be  an  optimal  choice  we  need  make  the  selection 


or 


p.(0/2,  <r) 


WP’ 


(1.30) 


P.(0/2.n) 


a+a* 


_?!  rYW  2iar(ep\ 

2a  +  2p\d»)  +«  +  /»W 

r1  \efi 1  r  «•  m\  a  ' 

+  *  [et'  +  PW  a  +  F\&) 


(l.?1) 


5.  Discussion  of  the  overall  reflection  and  transmission  coefficients 

With  an  explicit  p,  and  p0  at  our  disposal,  we  can  evaluate  the  fl(0,  <x)  and 
T(Q,o)  germane  to  our  elementary  trial  function.  For  reference,  their  complex 
variation  is  ploM«*d  as  a  function  of  the  discontinuity  angle  0  —  n,  in  Figure  2  for 
Poisson's  ratio  a  =  J.  The  present  magnitudes  |T|  are  somewhat  smaller  than  that 
given  by  earlier  first-order  calculations  (Kane  and  Spence  1963),  which  do  not 
simultaneously  yield  R  and  T.  Since  we  evaluate  both  R  and  T  together,  we  must, 
in  effect,  withdraw  some  energy  from  the  transmitted  field  to  allow  for  the  reflected 
wave.  Furthermore,  it  is  the  nature  of  the  variational  technique  to  underestimate 
the  subsidiary  diffraction  coefficients  p,  and  p„  since  it  only  yields  their  projection  in 
the  sub-space  spanned  by  the  trial  function. 

Although  the  analysis  is  certainly  valid  for  a  small  enough  discontinuity  in 
wedge  angle,  the  utility  of  the  procedure  cannot  be  established  until  there  is  some 
estimate  of  the  errors  committed.  A  feature  of  the  present  procedure  is  that  it 
suggests  a  natural  gauge  of  the  accuracy.  While  p,  and  p,  are  so  chosen  that  j]tj| 
and  \sj[  are  minimized,  both  sj  and  srr  are  non-zero  along  the  aperature  plane. 


R  MAG  ^  _T  MAG 
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These  residuals,  which  are  explicitly  known,  cannot  be  farther  reduced  without 
introducing  new  features  such  as  body  wave  contributions  into  the  analysis.  Since 
(I  -  |R|2  -  |T|J)  represents  that  fraction  of  energy  unaccounted  for,  we  can  estimate 
the  need  for  improving  the  calculations  by  examining  this  quantity. 

This  trror  estimate  is  a  eery  generous  one  because  only  part  of  it  implies  higher- 
order  corrections  to  R  and  7,'  the  remainder  representing  energy  which  is  accounted 
for  by  wave  to  body  wave  conversion.  The  data  of  Figure  2  show  that  if 

|0~n|  <  tO9, 

the  present  analysis  accounts  for  at  least  92  per  cent  of  the  energy,  and  therefore 
the  theory  should  not  require  further  improvement  within  this  range.  We  can  also 
compare  the  present  theory  with  experiment,  but  we  must  be  very  careful  if  we  do  so, 
because  there  are  fundamental  distinctions  betwccn'analysis  in  the  harmonic  domain 
and  pulse  measurements  (see  Appendix). 
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6.  Comment  on  the  formulation 

In  our  formulation,  we  have  used  as  boundary  conditions  the  demand  that  either 
the  normal  or  tangential  displacement  vanish  along  the  bisector  plane  depending 
upon  the  parity  of  tl»c  problem  under  consideration  Actually,  the  rigorous  bound¬ 
ary  conditions  along  the  plane  of  symmetry  are  somewhat  stronger  demands.  For 
example  in  the  symmetric  proM-tut  the  coirect  formulation  would  be  the  pair 
of  conditions  that  the  eompressioir.il  potential  be  an  even  function  about  the  plane 
of  symmetry  and  the  shear  potential  an  odd  function  about  tl.e  plane  of  symmetry, 
and  dee  versa  for  the  asymmetric  problem.  Then,  since  only  continuous  potentials 
arc  acceptable  as  solutions  within  the  wedge,  this  would  guarantee  that  there  would 
be  no  discontinuity  in  the  stresses  act  os  the  plane  of  symmetry.  For  an  approximate 
solution  however,  using  pairs  of  boundary  conditions  for  the  even  and  odd  problems 
is  unduly  complicated.  For  example,  consider  the  even  pioblcm;  any  trial  function 
under  consideration  would  have  shear  and  cotnpressionai  potentials  which  vary  at 
different  rates  along  the  bisector  plane.  The  boundary  condition  (1.18)  is  the  sum  of 
two  terms.  If  the  first  term  vanishes  then  (<£*  +  </>,)  !s  an  ei  c/t  function,  and  the  vanish¬ 
ing  of  the  second  term  implies  that  (t Pr  +  iIi,)  is  an  odd  functioa  As  a  result,  since 
we  minimize  (1.18)  in  the  mean  square  sense,  and  the  shear  and  congressional 
components  of  the  trial  function  vary  at  different  rates,  it  follows  that  minimizing 
the  sum  is  equivalent  to  minimizing  both  terms  in  the  sum  simultaneously. 


Fkx  3. — The  geometry  of  a  symmetrically  layered  wedge  which  can  support 
Love  waves.  In  the  figure,  it  is  assumed  that  the  height  H  is  such  that  only 
two  modes  propagate. 


Part  2.  Love  waves  on  an  clastic  wedge 

1,  Introduction 

A  layered  solid  can  support  surface  waves  which  are  horizontally  polarized 
shear  waves  trapped  in  the  superficial  layer.  Since  these  Love  waves,  as  they  are 
known,  have  no  compressions!  component,  it  is  not  necessary  to  introduce  potentials, 
and  it  is  possible  to  work  directly  with  one  scalar  function  w(x,  y),  the  z-component 
of  the  displacement  vector 


s  =  [0, 0,  wf.\\  y)) 


(2.1) 


Th*  theory  of  vurfjcr  wart  di'fruOion 


Within  the  £(  layer,  *v  satisfies  the  wave  equation 

(V1 + k>(.x,  y)  -  0,  k\  =  tu*Vt//»i.  (2-2) 

and  within  the  £a  substrate,  w  obeys 

(V1  +  k\Mx,  y)  =  0,  k\  =  o»Vj/^3-  (2-3) 

At  the  free  surface  of  £,.  the  normal  stress  must  vanish  and  this  will  be  true  provided 
that  the  normal  derivative 


vanishes  there.  ,  ,  .  .  , 

We  assume  the  £r£2  interface  to  be  welded  so  that  the  displacement  and 

normal  sires;  must  be  continuous  across  this  boundary.  Along  the  horizontal 
wedge  face,  those  conditions  will  be  satisfied  provided  that  the  Love  waves,  or 
jfrwavc s  of  amplitude  At  have  the  form 


A,  cos  J(k\  -  X})yeu“,  OZyZ-H  (2.5) 

A'y'  m  ^  A ,  cos#?  -XfH)  expJ{X>  -  A|)>  r4A'*,  y  Z  -  H,  (2-6) 

and  the  propagatioi.  '  ttants  Af  are  chosen  as  the  real  roots  or  the  period  equation 

(1,) 

If  the  shear  wave  is  to  be  trapped  in  the  layer,  or  ir  the  Love  wave  is  to  propagate 
we  need  Ik,!  £  U,  £  !fe,|-  For  any  thickness,  however  small,  there  is  at  least  one 
r^t  A? eorlefplndfng  an  acceptable  so.ution-the  fundamental  .5?, -wave.  AS 
the  acoustic  thickness  kjl  increases,  other  modes  can  propagate.  In  our  discussion, 
we  shall  assume  that  the  thickness  is  such  that  only  two  modes  propagate  the 
fundamental,  and  one  harmonic,  the  ^rwavc  with  a  propagation  constant 
The  analysis  proceeds  in  a  similar  fashion  if  an  arbitrary  number  of  modes  can 

propagate. 

2.  Formulation  of  the  boundary  value  problem 

We  assume  that  an  y,-wavc  is  incident  along  one  face  of  a  symmetrically 
layered  wedge.  At  the  discontinuity.  Tour  propagating  surface  waves  will  be  excited : 
a  reflected  and  transmitted  y,  wave  with  amplitude  coefficients  R„  and  7,, 
respectively,  and  reflected  and  transmitted  W 2  waves  whose  amplitudes  are  the 
conversion  coefficients  Rn  and  T,  2  respectively.  Our  task  will  be  to  determine  these 
ditTraction  coefficients  as  functions  of  the  wedge  angle  0,  the  layer  thickness  / 
and  the  elastic  constants  p„  p2.  ft,,  fh- 

By  the  same  argument  as  in  Part  1,  we  can  add  and  subtract  a  symmetric  J1,- 
excitation  on  the  light  wedge  face  which  leads  us  to  consider  a  pair  of  even  and  odd 
problems  in  a  bisected  wedge.  Since  Love  wave  d.ffraction  is  ?  scalar  problem,  the 
subsidiary  boundary  conditions  along  the  aperture  or  plane  ‘  ;>  mmetry  are  simp  y 

even  -  ^  =  0,  3  =  0/2— *  t2-^ 

r  d9 

odd  w  =  0,  9  =  0/2  — n  l2^) 


-  -  •  - - 
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for  the  even  and  odd  problems.  Tiie  trial  function  will  consist  of  an  incident 
,-wave,  and  retlected  and  y'J-wavcs,  with  unknown  amplitudes.  If  wc 

denote  the  subsidiary  reflection  and  conversion  coefficient  for  the  even  problem 
in  the  bisected  wedge  as  «•*,  and  ^2.  and  similarly  rf,  and  r*2  for  the  odd  problem, 
then  the  desired  major  coefficients  arc 

Kti  =  j(r*i  +  r" ,), 

(210) 

^12  *  i(r12  +  ru). 

(2.11) 

2*i l  “  K'*u-'1l). 

(212) 

^12  *°  i(r*  2  —  ^1*)- 

(213) 

As  in  Part  i,  wc  shali  determine  these  coefficients  by  a  variational  procedure  which 
ignores  body  wave  contributions. 


3.  Solution 

In  the  odd  problem  we  shall  choose  r?,  and  r*,2  so  that  the  residual  variation 
e(r)  along  the  aperture  plane 

e(r)  =  Srl+r*l¥*  +  S’llyi  9  «  0/2- n  (214) 

is  as  small  as  possible  in  the  mean  square  sense.  Using  the  same  definition  of 
scalar  product  as  in  Part  1,  we  have 

JeW!2  +  +  (215) 

and  this  expression  will  be  a  minimum  if.  and  only  if,  r*,  and  r{2  satisfy  the  normal 
equations 

&t,#l)+r*uWi,2>t)+r\i(£ri2’li  -  0,  (216) 

(ifl.J5fi)+r1lU'*,^-)  +  ^i(^?,^rJ)  =0.  (217) 

Equations  (216)  and  (217)  can  immediately  be  solved  for  r*u  and  r*u 


DErkir,.y,)  ur„?,t 


where 


1  (■2'„*\) 
DEr \wt,y2)  (y2,?2). 


DEr  = 


In  the  same  fashion,  for  the  even  problem  we  need  choose  r*,,  and  r*l2  as 


|fl££»  (l  e*'*  1  ggV\ 
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DET* 


(2.22) 


where 


DET* 


(\  02\  1  c!f\  ([  Ojn  [  Ojf\ 
j\r  09  '  r  09  )  \r  09  '  r  09  J 

(i  02\  \  oje_ A  A  ojrt  ay,\| 
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(2.23) 


With  this  knowledge,  the  reflection,  transmission  and  conversion  coefficients  arc 
given  by  (2.10)  to  (2,13). 


4.  Discussion  of  the  results 
(a)  Numerical  data 

We  have  used  the  preceding  formulas  to  calculate  the  diffraction  coefficients 
for  an  E.-layer  and  £,-substratc  for  which  k{/k2  =  1-297  and  =  2  159.  The 
!  phase  and  group  velocities  for  this  ease  have  been  given  by  Stoncly  and  are  available 
in  a  standard  reference  (p.  213,  Ewing  &  others  1957).  Figures  4  to  7  illustrate  the 
variation  of  the  magnitude  of  the  diffraction  coefficients  which  are  even  functions 
of  the  discontinuity  angle  ©-«.  The  curves  arc  indexed  by  four  values  or  the 
dimensionless  parameter  A,H,  namely  5, 6, 7.8;  ir  A,H  =  5.  then  the  second  mode  is 
just  above  cut-off  and  if  x.tf  =  8,  the  third  mode  is  just  below  cut-off.  It  is  very 
interesting  to  note  that  the  conversion  coefficient  Tl2  exceeds  the  reflection  coeffi¬ 
cient  R, That  is,  there  is  a  tendency  for  the  energy  to  continue  to  propagate  in  the 
same  direction  even  if  it  necessitates  a  transfer  of  model  characteristics. 

(b)  Interpretation 

If  we  compare  any  two  waves  of  identical  characteristics,  then  their  relative 
energy  is  proportional  to  the  absolute  square  or  any  corresponding  amplitude.  On 
the  other  hand,  before  we  can  compare  the  energy  in  the  fundamental  2  ,-  vase  to 
that  of  an  y2-wavc,  its  first  harmonic,  we  n:ed  make  some  further  calculations. 
With  no  loss  in  generality,  let  us  specialise  our  uiseussion  to  the  horizontal  wedge 
face  for  which  the  y,-w  ,.cs  arc  given  explicitly  by  {2.5)  and  (16),  and  evaluate  t  e 
scalar  product  along  the  wavefront  y  =  0.  Thus  \A,\2  jjy,i!‘  represents  the  mean 
square  energy  flux  transported  by  an  y.-wave  of  amplitude  A>.  If  we  denote  the 
group  velocity  of  an  yrwave  as  A,  it  follows  that  the  ratio 

Kl’A.jy.E* 


l^pA^y-ii2 


(2.24) 


compares  the  power  flow  of  an  /l.y.-wavc  to  that  of  an  42y2-wavc.  In  particular, 
a  mode  near  cut-off  behaves  like  an  unbounded  plane  wave  m  the  Ej-mcdium; 
hence  such  a  wave  can  carry  large  amounts  of  power  even  if  its  amplitude  is  decep¬ 
tively  small.  As  a  result,  if  we  are  to  discuss  power  transfer,  we  should  renormalize 
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Fig  4.— The  magnitude  of  the  Love  wave  diffraction  coefficients  for  a,H  -  5 ; 
they  are  even  functions  of  the  discontinuity  aigle  0-n.  In  this  case 
Ik,/*,  -  1-297 ■  2159;  the  normalization  value  N  -  2-0. 


the  amplitude  of  the  conversion  coefficients 

Ri i  ~  ‘  ii  -  NTtl 

where 

_  Ai  i *A 

~  a,  Fsr,| 


(2-25) 

(2.26) 


so  that  |R'j2|j  and  |Tf  j|J  arc  proportional  to  the  power  transferred  by  the  diffraction 
of  a  if, -wave  of  unit  amplitude  at  a  .edge  discontinuity.  The  appropriate  values 
of  N  for  the  previous  numerical  example  aYc  cited  in  the  captions  of  Figures  4  to  7. 
in  a  fashion  similar  to  the  error  analysis  of  Part  1,  the  function 


(i-jR„  MW-F.ilHW) 


represents  the  amount  of  ambiguous  energy.  These  values  arc  much  more  satis¬ 
factory  in  the  present  analysis  than  in  Part  1.  For  example,  when  ktH  *=  7  and 
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Abso  uie  Discontinuity  Angle  |0-u|  Absolute  Discontinuity  Angle  J©-*!  — 

magnitude  of  the  Love  wave  difTr.-ttion  coeff  ents  for  Ft&  6.— The  magnitude  of  the  Love  wave  diffraction  coefficients  for 

are  even  functions  of  tne  discontinuity  angle  0-r.  In  this  ltH  -  7;  they  are  even  functions  of  the  discontinuity  angle  O-ic.  In  this 

.  1.707  »  71  SO ■  the  normalization  value  N  M  case  kJk.  »  1297,  ujit.  •=  2150.  the  normalization  value  Af  «  I  I 
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Fi&  7.— The  magnitude  of  the  Love  wave  difTraction  coefficients  for 
A,H  8;  they  arc  even  functions  of  the  discontinuity  angle  0—  it  In  this 
case  fcj/fc,  -  1.297,  nl!nt  *»  2.159;  the  normalization  value  N  -  1.0. 


the  discontinuity  angle  |0-nj  is  as  large  as  30\  the  present  analysis  accounts  for 
at  least  95  per  cent  of  the  scattered  energy.  This  improvement  is  due  in  large  part 
to  a  use  of  a  more  flexible  trial  function.  In  Part  2.  'wo  trial  coefficients  have  been 
available  for  either  symmetry. 
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Appendix 

Love  wave  diffraction  coeflkients  would  be  very  difficult  to  measure  in  the 
laboratory,  but  the  techniques  of  two-dimensional  model  seismology  offer  a  means 
of  determining  Rayleigh  wave  reflection  and  transmission  coefficients  with  an 
accuracy  ol  about  10-20  per  cent.  The  present  thcc*y  and  experiment  agree  if 
0  -  n,  but  outside  this  range  there  arc  experimental  features  which  arc  not  dupli¬ 
cated  by  the  results  of  the  present  elementary  variational  procedure.  The  analysis 
in  the  harmonic  domain  could  be  refined  by  employing  various  devices  to  reduce 
the  amount  of  unexplained  energy.  Such  calculations  would  probably  require 
substantial  effort,  and  the  idealized  formulation  of  lie  present  problem  should  be 
reviewed  if  the  labour  is  to  have  relevance  to  pulse  rncasur  ments. 

There  arc  major  distinctions  between  analysis  in  the  harmonic  and  time  domain. 
For  example,  whereas  a  harmonic  Rayleigh  wave  is  a  uniquely  defined  entity, 
Friedlander  (194S)  has  pointed  out  that  a  Piylcigh  pulse  can  assume  a  variety  of 
waveforms.  Furthermore,  any  Rayleigh  pulse  cannot  have  a  sharply  defined  wave- 
front,  and  theoretically  must  give  infinite  advance  notice  of  its  arrival  :nless  it 
merges  continuously  with  a  precursor,  typically  the  shear  pseudo-surface  wave 
(Cagniard  1939).  Although  the  amplitude  o'  this  shear  wave  decays  with  distance, 
its  integrated  flux  remains  constant.  It  is  difficult  to  separate  the  far-ficld  effects  due 
to  the  arrival  of  a  Rayleigh  pulse  and  its  shear  companion  at  the  second  wedge  face. 
In  other  words,  in  addition  to  Raylcigh/Ravlcigh  interactions,  there  will  be  some 
shcar/Raylcigh  conversions.  What  contributions  might  this  shear  wave  introduce? 
Whereas  we  cannot  give  a  rigorous  answer  to  this  question,  we  can  however  make 
a  rough,  but  simple,  estimate. 

We  first  note  that  if  the  wedge  angle  0  is  n  oi  a/2,  then  we  would  expect  little 
or  no  shear/Rayleigh  conversion.  In  ti  e  first  ease,  there  is  a  discontinuity,  and  the 
second  case  corresponds  to  a  geometry  for  which  the  shear  wave  is  essentially 
normal  to  the  second  wedge  face,  and  we  know  that  for  normal  incidence,  a  shear 
wave  is  reflected  as  a  shear  wave.  Then,  from  Equation  (1.6),  we  note  that  we  can 
easily  add  some  additional  shear  potential  to  the  original  excitation  by  incrementing 
the  Rayleigh  wave's  shear  coefficient  T  by  an  additional  contribution  /(©)  depend¬ 
ing  on  the  wedge  angle 


r  =  r[i  +/(©)]. 
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In  this  stud.  we  describe  proce Jure*  for  the  numerical  solution  of  the  sccoml-ordcr  linear  differential 
equation  wliieli  have  either  continuous  or  discontinuous  coefficient*.  Our  motivation  is  a  well-known 
technique  in  the  theory  of  inhomogeneous  transmission  lines  r  the  treatment  of  a  continuously  varying 
line  by  considering  it  to  be  comprised  of  various  sections  of  uniform  lines.  Although  this  procedure  is 
very  suggestive  physically,  there  are  difficulties  with  its  applications  to  sccond-oidcr  equations  whicii 
do  not  describe  wave  propagation.  First,  the  language  of  the  circuit  engineer  is  auch  that  a  pair  of 
first-order  equations  describing  some  analogs  of  the  complex  quantities  voltage  and  current  seem  to 
be  required.  Whereas  these  equations  appear  naturally  in  transmission  lino  theory,  we  show  that  it  is 
an  unnecessary  burden  to  find  their  counterparts  when  the  problem  is  hut  to  solve  a  second-order 
equation.  The  second  objection  to  this  approach  is  that  it  is  not  apparent  tliat  a  piecewise  constant 
partition  of  the  coefficients  in  a  differential  equation  will  yield  the  rigorous  solution  if  the  subdivision 
is  canied  out  to  an  arbitrary  degree.  Indeed,  we  show  that  the  limit  of  the  quantitation  scheme  can 
not  yield  the  rigorous  solution.  On  the  other  hand,  wc  are  led  to  a  well-defined  technique  which 
generates  the  solution  by  a  method  suggested  by  the  procedures  used  in  the  discrete  case.  This  enables 
the  second-order  equation  to  be  solved  rigorously  by  iteration  with  no  complications,  and  in  a  form 
ideally  suited  for  computer  programs. 


L  INTRODUCTION 
A.  Basic  Equations 

HE  class  of  problems,  mathematical  and  phys¬ 
ical,  which  require  an  explicit  knowledge  of  tho 
solutions  of  the  linear  second-order  differential  equa¬ 
tion 

d'U(t)/dt'  +  A  (t)dU(t)/dt  +  B(t)  U(t)  -  0  (1) 

is  legion.  Only  rarely  can  exact  solutions  of  this 
equation  be  found  in  closed  form.  While  local  approx¬ 
imations  are  easily  constructed,  the  task  of  finding 
global  representations  offers  a  major  problem.  This 
difficulty  is  particularly  vexing  inasmuch  as  the 
first-order  equation, 

dV/dt  -  P(l)V(l)  -  0  (2) 

possesses  an  immediate  integral 

V(t )  -  exp  {/'  P(r)  dr}  .  (3) 

but  no  similar  representation  for  Eq.  (1)  has  been 
found.  The  solution  (3)  has  the  interpretation  of 
defining  the  present  value  of  K(f)  in  terms  of  the 
history  of  its  logarithmic  derivative  P(f)  = 
(d/df)(In  V).  If  wc  specify  an  initial  conditiou  on 
the  possible  pair  of  solutions  of  (1)  at  say  t  «  f„ 
then  Eq.  (1)  describes  a  unique  function  with  that 
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property.  From  this  initial  condition  and  a  knowl¬ 
edge  of  Eq.  (1),  the  behavior  of  the  second-order 
!/(/)  should  likewise  be  discernible  in  terms  of  the 
prior  trajectory  of  its  logarithmic  derivative:  the 
topic  of  the  present  analysis. 

Following  a  change  of  unknown 

u(0  -  U(l)  exp  ||  J'  A(t)  dr  j  (4) 

and  an  introduction  of 

*’(f)  -  BO  -  idA«)/dt~  H’(0,  (5) 

Eq.  (1)  can  be  recast  into  the  canonical  form 

d*u/dx’  +  Jt*(i)u(i)  -  0  (6) 

in  which  we  replace  f  by  x  to  eliminate  confusion 
between  the  forms  (l)  and  (6).  At  some  point  x  —  z%, 
wc  may  suppose  that  we  have  a  given  constraint 

5  b  u  “  "  0*  (*  -  *•).  (7) 

which  singles  out  one  solution;  another  linearly  in¬ 
dependent  solution  can  be  specified  by  choosing  a 
different  value  for  the  parameter  g9. 

The  form  of  the  solution  (3)  to  the  first-order 
equation  (2),  suggests  that  it  might  be  more  profit¬ 
able  to  concentrate  our  attention  on 

/  r(f)df  -  lnu(x),  u(x)  -  exp  /  r({)  d$,  (8) 

•»» 
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that  is,  the  logarithmic  dcri  'ative 

r(i)  -  d4~ . 


(9) 


rather  tiian  on  u(i).  Indeed,  from  the  point  of  view 
of  a  mathematical  physicist  interested  in  eigenvalue 
problems,  an  explicit  knowledge  of  the  logarithmic 
derivative  provides  more  insight  than  would  form¬ 
ulas  for  two  distinct  solutions  of  (l),  therefore  wc 
shall  pay  special  attention  to  the  logarithmic  deriv¬ 
ative.  The  continuity  of  P  is  assured  if  wc  restrict 
our  attention  to  continuous  solutions  u(z)  possessing 
a  continuous  first  derivative  u‘(x).  If  these  functions 
are  bounded,  then  rvc  note  that  the  poles  and  zeros 
of  r(x)  correspond  to  the  zeros  and  maxima  re¬ 
spectively  of  u(i).  In  terms  of  r(x),  the  equations 
(6)  and  (7)  become 


dr/dx  ~  -(it*  +  T*),  x„  <  x  <  x. 


r  “  ga, 


X  ■=  x» 


(10) 

(U) 


wherein  we  have  adopted  the  convention  that  the 
initial  condition  is  given  at  x  «*  xa,  and  wc  wish 
to  propagate  the  solution  to  the  right  of  x*  up  to 
Borne  arbitrary  point  x„.  Formula  (10)  describes  T(x) 
in  terms  of  a  Ricatti  equation,  and  even  though 
this  expression  is  nonlinear,  it  can  easily  be  solved 
by  iterative  techniques  if  fc’(x)  is  negative.  On  the 
other  hand,  if  fc*(x)  is  positive,  then  T(x)  behaves 
like  an  ill-tempered  function  with  numerous  sing¬ 
ularities.  Even  high-speed  computers  have  difficulty 
in  coping  with  such  a  wildly  fluctuating  function, 
and  numerical  work  under  such  circumstances  oflc-s 
major  difficulties.  One  conclusion  of  the  present  work 
is  the  need  to  isolate  and  describe  explicitly  the 
bad  manners  of  T(x)  so  that  if  numerical  methods 
must  be  used,  they  should  treat  slowly-varying 
functions. 

B.  Summary  of  Results 

Equation  (6)  can  be  thought  of  as  a  reduced  wave 
equation:  subject  to  severe  limitations  on  the  be¬ 
havior  of  Jt’(x),  approximate  solutions  can  be  found 
by  the  phase-integral'  method.  This  technique  yields 
reasonably  simple  solutions  of,  say,  exponential  type. 
However,  the  character  of  the  solution  must  change 
in  the  vicinity  of  a  turning  point  of  fc’(x):  in  such 
a  region,  the  solution  must  be  oscillatory  rather  than 
exponential  which  means  that  the  phase  integral 
nee'1 3  severe  correction.  Loosely  speaking,  the  major 

1  We  prefer  this  adjective  to  any  acronym  based  upon 
tome  perturbation  of  the  initials  of  .Uouvtlle,  Green,  Carlini, 
Stokes,  Horn,  Rayleigh,  BirkhofT.  Jeffreys,  Langcr,  Wentrel, 
Kramers,  Briilouin,  tl  al. 


disadvantage  of  the  phase-integral  method  is  that 
it  involves  an  argument  of  the  form  /*  k(i)  in 
trigonometric  or  exponential  functions.  As  a  result, 
wl  en  fc’(x)  changes  sign,  the  argument  docs  not 
change  character  correctly  because  of  the  integral’s 
inertia  in  storing  the  past  values  of  k(x).  That  is, 
the  phase-integral  solution  must  involve  complex 
quantities  in  the  presence  of  turning  points  even 
if  the  original  differential  equation’s  codhcicnts  are 
real.  We  know  that  such  equations  have  real  solu¬ 
tions,  but  the  phase  integral’s  hysteresis  prevents 
its  providing  a  valid  representation  of  them.  For 
this  reason,  we  avoid  using  formulas  which  depend 
upon  integrals  of  the  form  /*  Jfc(()  </{. 

One  means  of  avoiding  turning-point  difficulties 
is  to  consider  k7(x)  as  being  piecewise  constant 
i*(x)  *  k]  over  sufficiently  small  intervals  x,_,  < 
x  <  x4.  Within  each  quantum  cell  a  pair  of  trial 
solutions  can  be  inserted:  trigonometric  or  expo¬ 
nential  depending  upon  the  sigir  of  JkJ.  By  this 
process,  the  task  becomes  an  algebraic  one  of  de¬ 
termining  the  correct  coeflicicnfs:  of  the  trial  function 
in  each  ccl.  so  that  it  joins  smoothly  with  its  neigh¬ 
bors.  The  difficulties  with  this  procedure  are  twofold. 
First,  while  the  procedure  appears  to  converge  to 
a  solution  of  the  basis  of  physical  arguments,  the 
mathematical  justification  for  this  procedure  has 
neither  been  simple  nor  convincing,  since  the  limiting 
process  involves  a  sequence  of  step-function  approx¬ 
imations.  Second,  for  fine-grained  quantization,  the 
large  number  of  algebraic  equations  to  be  solved 
can  be  overwhelming.  In  diffraction  theory,  an  ap¬ 
proximate  means  of  surmounting  this  hurdle  is  the 
heuristic  device  of  neglecting  certain  contributions 
within  each  cell  that  can  be  identified  with  multiple 
reflections  of  a  wave  at  a  boundary  which  has  a 
very  small  reflection  coefficient.’ Indeed,  the  phase 
integral  is  but  the  leading  term  in  such  an  expan¬ 
sion.'  In  Sec.  II,  wc  consider  this  discretization 
process  .nd  show  how  the  algebraic  difficulties  can 
be  sidestepped  by  introducing  the  proper  variables 
and  format,  of  the  trial  solution  without  making  any 
assumption  concerning  the  behavior  of  the  field 
within  each  cell.  This  yields  an  expression  for  I\x) 
in  terms  of  a  product  of  n  2  X  2  unimodular  matrices 
for  an  nth-order  quantization.  In  the  limit,  as  n  — < »  <*> , 
each  matrix  operator  becomes  an  infinitesimal  trans- 
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formation;  \vc  analyze  this  differential  change  in  E 
and  find  a  geometric  characterization  of  the  tra¬ 
jectory  of  the  logarithmic  derivative.  We  then  find 
that  the  quantization  scheme  cannot  yield  the 
rigorous  solution  if  the  coefficients  in  the  original 
differential  equation  arc  analytic.  Nonetheless,  the 
procedures  suggest  a  transformation  which  cnahlcs 
the  rigorous  solution  to  be  found  by  a  related  tech¬ 
nique.  This  analysis  will  be  found  in  See.  Ill  together 
with  an  illustrative  numerical  example. 

II.  DISCRETE  FORMALISM  AND  SOLUTION 

A.  Introduction 

\\  c  assume  that  gt  and  4’(x)  are  real  [but  not  4(x) j 
so  that  we  can  confine  our  attention  to  real  solutions 
of  (6).  Of  course,  once  the  real  solutions  of  (G)  arc 
known,  complex  solutions  can  easily  be  constructed 
by  the  superposition  principle. 

Let  [4,|"  be  any  monotone  sequence  of  (m  +  1) 
points  which  divides  any  portion  of  the  x  axis,  say, 
the  interval 

<*<*.  <+•  (12) 

into  subintervals 

Ax,  -  a,  -  a,.,,  (a,  -  X„  ai4t  >  a„  a„  -  xj.  (13) 

If  this  quantization  is  so  fine  that  within  the  »th 
interval  Ax,  we  can  replace  fc’(x)  by  fc’  to  any  pre¬ 
assigned  degree  of  accuracy,  then  the  equation 

(fut/dx1  +  k*  u,  -  0,  <  *  <  a,  (W) 


and  so  ,1 ,  and  T ,  can  be  identified  as  the  amplitude 
and  logarithmic  derivative  of  the  solution  at  the 
left  endpoint.  Most  importan'  however,  is  that  the 
form  of  (15)  is  such  that  both  the  boundary  condi¬ 
tions  (17)  are  decoupled,  and  it  becomes  possible 
to  solve  for  Ai  and  r,  separately:  a  considerable 
reduction  of  computational  effort. 

B.  Solution 

The  equations  (G)  and  (7)  arc  homogeneous,  and 
so  we  can  normalize  A a  to  any  convenient  value, 
but  we  must  choose  T0  =  gt  if  the  boundary  condi¬ 
tion  (11)  is  to  be  satisfied.  Even  though  4,  is  a  jump 
function,  this  initial  data  can  be  continued  as  a 
smooth  function  for  all  x,  provided  we  choose 

T  -  r<  k,Ax,  -  k,  sin  4.  Ax, 

T >K  sin  4, Ax,  -f-  cos  k.  Ax,  '  ' 

and 

Aitt  *=  j4,[r,A*(  ‘  sin  I, Ax,  +  cos  4, Ax,].  (19) 

In  this  manner  we  obtain  a  function  w  hich  satisfies 
(6)  everywhere  except  that  4*(x)  must  be  replaced 
by  some  local  average  value  It J.  Furthermore,  the 
iterations  required  to  extend  u,(x)  are  very  cosily 
performed  since  (18)  has  the  form  of  a  Mobius 
transformation  whose  coefficient  mat.  ix 

M,  -  cosI.Ax,  -4.sinI.Ax,  ^ 
fiV  si»  k,  cos  4,  Ax, 


has  the  general  solution 

v,  •»  .4,[r,I-r‘sir  4,(x  —  4..,) 

+  cos  k,(x  -  3,_t)J,  (15) 

wherein  the  discrete  parameters  A ,  and  T,  remain 
to  be  specified.  From  a  theoretical  point  of  view, 
any  two-parameter  combination  of  sin  ifc,i  and 
cos  4,x  could  serve  as  well  as  the  particular  form 
selected  in  Eq.  (15).  However,  from  a  practical  view¬ 
point,  major  advantages  arc  introduced  by  the  selec¬ 
tion  (15).  First,  no  assumption  has  been  made 
concerning  the  character  of  u,  in  the  ith  interval: 
oscillatory’  or  damped.  For,  should  k\  be  negative, 
or  4,  «  it,  be  imaginary,  then 


u,  *  .4,[—  r.*;’  sinh  *,(x  —  4,_,) 

+  cosh*,(x  —  4,_,)j  (ifi) 

w~ll  continue  to  be  real  provided  only  that  .1,  and 
T,  arc  chosen  real.  Note  that  at  x  -  4,.,  we  have 


haa  its  determinant  equal  to  unity.  As  a  result, 
the  coefficient  matrix  T„  required  to  transform  I\ 
into  some  final  value  r.  at  the  arbitrary  point 
the  product  of  m  factors 

T„  •  ...  -M,.  ...  (21) 

and  is  likewise  unimodular. 


C.  Check 

It  is  now  but  a  simple  calculation  to  show  that 
if  m  — »  w,  and  Ax,  — *  0  for  any'  sequence  (4,)", 
then  the  process  yields  a  function  which  converges 
to  the  solution  of  an  equation  consistent  with  as¬ 
sumed  behavior  of  k\x).  Since  T(x)  defines  u(x) 
uniquely  apart  from  a  multiplicative  constant,  it 
will  be  sufficient  .o  show  that  the  continuous  limit 
of  the  r,’s  defined  (18)  satisfies  the  Ricatti  equation 
(10)  on  an}'  interval  Ax,.  We  calculate  that 


_-4-,  sin  I, A.r,  -  (r?4~')  sin  4.  At. 
^J-.[(r,ir')  sin  4. Ax,  +  coe  4, Ax.)  ’ 
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and  if  we  pass  to  the  limit  Ax,  — *  0  we  obtain 

Iim  r<,,A— -  V;  **  -a*  +  r’),  (23) 

a., -a  « 

which  is  precisely  Eq.  (10).  As  a  double  check,  we 
can  repeat  the  limiting  process  for  Eq.  (19), 

A,.,  -  A,  __  A.tfr.t;1)  sin  A|Ax,  +  cos  t.Ax,  -  I) 
Ax.  Ax, 

(24) 


and  derive  the  defining  equation  for  T(r)  in  terms 
of  m(j) 


A , . ,  —  A ,  du 
**  dz 


Iim  — 


-  «r(x). 


(25) 


if  we  note  that  in  the  limit,  the  amplitude  A, 
becomes  u(x,). 


Since  we  have  always  assumed  that  dk/dx  *»  0 
locally,  what  we  have  just  shown  is  that  to  the  extent 
that  the  derivative  k'  can  be  neglected  on  any  quantum 
cell  the  scheme  converges  to  a  solution  of  (C). 


age  and  current,  or  E  and  H  are  out  of  phase,  and 
the  complex  entries  in  the  transmission-line  matrix 
emphasize  tins  feature.  In  our  analysis  however,  we 
have  avoided  introducing  any  analogs  of  the  first- 
order  transmission-line  equations  since  we  feel  that 
their  use  i.i  artificial  and  contrived  for  applications 
other  than  those  described  by  Maxwell’s  equations. 
Furthermore,  unlike  the  lh,‘s.  their  geometric  in¬ 
terpretation  requires  a  complicated  hyperbolic  model 
rather  than  the  more  familiar  Euclidian  plane. 

We  emphasize  that  Eq.  (20)  represents  a  real 
transformation:  in  more  genera!  terms,  it  represents 
a  mapping  of  the  complex  T  plane  upon  itself  such 
that  the  real  axis  is  left  invariant.  Since  Mdbius 
transformations  describe  elementary  geometric  de¬ 
formations  of  the  complex  plane,  this  pro. Ides  us 
with  a  clue  to  the  proper  characterization  of  M,. 
By  direct  multiplication,  it  is  easily  checked  that 
each  M,  can  be  re-expressed  as  the  product  of  three 
terms 

M,  -  Cr'R.C,,  (26a) 


D.  The  Infinitesimal  Transformation 

The  preceding  analysis  furnishes  a  simple  illustra¬ 
tion  of  the  algebraic  structure  of  a  particular  Ricatti 
equation:  one  example  of  the  general  theory  of  such 
equations*  w-hich  shows  that  its  solution  can  be  char¬ 
acterized  by  matrix  operators.  In  the  limit,  each 
matrix  M,  in  the  product  (21)  approaches  an  infini¬ 
tesimal  transformation,  and  the  total  transformation 
T.  giving  T(z.)  in  terms  of  r(x0)  can  be  expressed 
either  as  a  matrizant7  or  a  product  integration*:  the 
finite  result  of  compounding  a  sequence  of  infinites¬ 
imal  transformations.  We  can  think  of  this  operation 
as  defining  a  curve  traced  out  by  the  repeated 
reaction  of  the  original  point  T(x0)  to  each  infinites¬ 
imal  product  M,.  If  we  can  find  some  analytic  means 
of  specifying  this  trajectory,  then  we  might  be  spared 
the  need  of  calculating  matrizants,  or  product 
integrals  to  describe  the  limit  of  the  iterative  process. 
For  this  purpose,  we  consider  the  geometry  of  the 
transformation  implied  by  any  matrix  M(. 

No  doubt  sonic  readers  have  been  struck  by  the 
resemblance  of  M(  to  the  transmission-line  matrices 
of  network  theory.  However,  we  caution  that  this 
resemblance  is  deceiving:  M<  is  a  real  transformation 
unlike  the  transmission-line  matrix  which  is  a  complex * 
one.  In  the  context  of  transmission-line  theory,  volt- 

•  R.  RedhelTcr,  i.  Rati.  Mech.  Anal.  5,  835  (195G). 

’  R.  A.  Frazer,  W.  J.  Duncan,  and  A.  R.  Collar,  Elt- 
mtntnry  Mafrirtt  (Cambridge  Univenity  Prm,  New  York, 
1959). 

•  G.  Birkhoff,  J.  Math  A  Phyaiw,  IS,  104  (1937). 


where 

V*  0 

C,  -  u 

.  0  it:', 

cos  k,  Ax4  -sin  k.Ax, 

K,  ■* 

sin  k, At4  cosk.Ax,. 

Each  fuctor  lias  a  simple  geometric  interpretation. 
The  first,  C,  represents  a  uniform  contraction  of 
the  r'-plauc  by  a  factor  1/Jt(;  on  the  other  hand, 
its  inverse  C7*  represents  a  uniform  dilation  by  the 
inverse  factor  k,.  Between  these  two  operations,  we 
need  perforin  R,  which  represent?  a  rotation  of  the 
stcrcographic  projection  of  the  T'-plnnc  on  the 
Riemann  sphere. 

The  sequence  of  operations  is  illustrated  in  Fig.  1 
which  only  shows  the  plane  intersecting  the  real 
T-axis  and  the  north  pole  N  o f  the  Riemann  sphere. 
First,  we  lay  of!  the  initial  point  T«  on  the  T-axis 


Flu.  I.  Tbc  geometry  ot  the  matrix  traiuformatiuu  M, 


(26b) 

(26c) 
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which  has  been  drawn  obliquely  for  convenience. 
The  r'-  axis  is  drawn  horizontally  and  intersects  the 
F-axis  at  an  angle  cos'1  (1  /k,).  If  we  project 
the  T-axis  onto  the  F'-axis,  then  all  lengths  on  this 
new  axis  will  be  uniformly  contracted.  We  continue 
by  finding  the  stcrographic  projection  of  this  image 
point,  and  then  we  rotate  the  Ricinnnn  sphere  by 
a  negative  (clockwise)  angle  -6,  ■=  tan"1  (fc.Ar.)  if 
Jt.Aar,  is  positive.  We  retrace  our  steps  to  the  original 
T-axis,  and  it  is  a  simple  matter  to  show  that  the 
aforementioned  sequence  is  mimicked  by  the  matrix 
product  M.  =»  C-R,C(.  The  proof  of  this  assertion 
is  left  to  tho  reader,  and  is  easily  reproduced  once 
it  is  recalled  that  the  angle  subtended  by  an  arc 
ou  a  circle  from  the  center  is  twice  the  angle  sub¬ 
tended  by  that  arc  from  the  north  pole. 

However,  the  geometric  construction  also  shows 
that  any  transformation  AT,,  can  be  obtained  by 
differentiating  the  function  I\  «  kt  tan  (—k,x), 
that  is 

AT,  =  (d/dx)[ki  tan  (-fc,x)l,.„Ax,.  (27) 

If  we  recognize  that  this  transformation  must  evolve 
continuously  from  somo  initial  value  F,  -  g,  at 
x  «  X{,  we  thus  find  the  general  integre1 

T(x)  -  k,  tan  [k,(y<  -  x)]  (28) 

representing  the  locus  of  r(x)  on  the  tth  interval 
where  the  constant  y,  is  fixed  by  the  initial  condition 
at  x  «  r 

7.  -  x,  -  Jtr’  tan"  (0i/k.).  (29) 

In  the  limit,  as  Ax*  — *  0  it  would  seem  that  k,  could 
be  replaced  by  A:(x)  and  g,  by  its  initial  value  g„ 
since  g,  —  gt.  x  — *  0.  In  other  words,  we  expect 
that  the  net  transformation  could  be  described  as 

T(x)  =  k(x)  tan  [(y  -  x)*(x)l 

7  “  x.  +  (i/*(idj  tar'  :  ig«A(x,)),  (30) 

where  g0  is  the  initial  value  of  F  at  x  ■=  xt.  However, 
the  promise  of  this  anticipation  is  broken  for  if  we 
substitute  (30)  into  the  Ricatti  equation  we  find. 

—  -  -(**  +  T1)  +  (dk/dx)  {tan  [(y  -  x)k\ 

+  k(y  -  x)  see’  [(y  -  x)t]|.  (31) 

In  other  words,  (30)  satisfies  the  Ricatti  equation 
only  to  the  extent  that  the  first  derivative  of  k(z) 
can  be  ncg’°cted.  Where  is  the  fallacy  in  the  preced¬ 
ing  argument  that  can  explain  our  disappointment, 
and  perhaps  lead  us  to  a  correction? 
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The  resolution  of  this  conundrum  depends  upon 
a  subtle  point  in  the  character  of  the  analysis.  The 
basic  assumption  was  that  a  continuous  k'(x)  could 
bo  approximated  arbitrarily  well  by  a  staircase 
approximation.  Whereas  any  such  approximation  is 
discontinuous,  it  does  approach  ft  continuous  limit 
as  the  treads  and  risers  get  finer  and  finer.  Of  the 
eventual  continuity,  there  is  no  doubt,  but  it  is  ft 
peculiar  type  of  continuity.  Think  for  example  of  a 
staircase  approximation  to  k’(x)  “ion  the  interval 
0  <  x  <  1.  Each  riser  represents  the  jump  from  one 
tread  to  the  next,  and  the  height  of  any  one  riser 
goes  to  zero  as  the  approximation  gets  better  How¬ 
ever,  no  matter  how  fine  an  approximation  we  care 
to  make,  the  sum  of  all  the  users,  or  the  total  jump 
xoill  always  be  constant.  In  precise  terminology,  this 
is  the  difference  between  continuity  and  absolute 
continuity.  The  fallacy  is  now  evident:  while  y,  — 
Yi-i  -♦  0,  we  can  not  conclude  that  y,  will  be  un¬ 
changed  as  we  pass  through  an  unbounded  number 
of  quantum  cells.  That  is,  as  we  proceed  from  the 
tth  cell  to  the  (t  -f  l)th  one,  we  can  not  neglect 
the  differential  jump  y4»,  —  y,  at  the  end  of  the 
ith  interval.  In  other  words,  we  can  think  of  Eq.  (30) 
as  a  representation  of  the  continuous  or  principal 
part  of  the  transformation  (21)  which  overlooks  the 
denumerable  collection  of  minute  jumps. 

We  can  now  see  that  the  quantization  scheme 
can  not  yield  the  rigorous  solution  even  if  carried 
out  to  the  continuous  limit,  for  the  result  of  a 
staircase  approximation  to  Jfc*(x)  can  not  be  ab¬ 
solutely  contin','"'s.  From  the  differential  equation 
,r\  ...at  the  ratio  u"/u  must  be  equal  to 

—fc’tx),  hence  m  the  limit,  u"/u  can  not  be  ab¬ 
solutely  continuous  no  matter  how  smooth  the 
original  k2{x)  might  have  been.  In  general,  we  know 
that  whenever  i’(x)  is  analytic,  the  differential  equa¬ 
tion  (6)  must  have  an  analytic  solution.  As  a  result, 
for  such  equations  the  fact  that  the  limiting  u”/u 
is  not  absolutely  continuous  contradicts  the  required 
analyticity  of  the  solution.  As  a  result,  we  conclude 
that  the  limit  of  the  quantization  process  in  general 
will  not  converge  to  the  rigorous  solution.  That  the 
limit  of  the  iterative  process  might  be  a  weak  solu¬ 
tion  in  terms  of  some  norm  remains  an  open  question. 
On  physical  grounds,  there  can  be  no  doubt  that 
the  discretization  approach  is  a  useful  procedure 
when  k*(x)  can  be  considered  to  be  well  approximated 
by  a  piecewise  constant  function.  Ferhaps  if  we 
retain  information  concerning  the  derivatives  of 
fc’(x),  we  can  construct  a  rigorous  solution?  This 
question  is  answered  affirmatively  in  the  next 
section. 


-  if? 
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HI.  THE  CONTINUOUS  TRANSFORMATION 

In  the  preceding  paragraphs,  wc  have  discussed 
the  hazards  that  obstruct  an  extension  of  the  finite 
quantization  theory  to  the  continuous  ease.  In  such 
&  situation,  it  is  always  well  to  pause,  und  reflect 
upon  the  possiblity  that  the  infinitesimal  case  might 
better  bo  handled  by  other  methods.  Basically,  wc 
arc  considering  the  second-order  equation  as  an 
initial-value  problem.  That  is,  wc  wish  some  means 
of  knowing  the  direction  the  solution  or  its  log¬ 
arithmic  derivative  will  take  once  an  initial  value 
is  specified  This  question  is  of  course  immediately 
answered  by  the  Kicatti  equation  (10)  for  it  describes 
the  derivative  cxplicity  in  terms  of  ‘lie  function. 
However,  wc  know  that  when  *s( x)  is  positive, 
r(x)  behaves  like  the  tangent  function  and  has 
numerous  poles  and  zeros.  Any  numerical  method 
which  attempts  to  pursue  the  path  of  such  a  spirited 
function  will  almost  of  necessity  be  doomed  from 
the  outset.  On  the  other  hand,  the  poles  and  zeros 
of  the  logarithmic  derivative  have  no  more  signif¬ 
icance  than  locating  the  maxima  and  zeros  of  the 
desired  solution  of  the  differential  equation  and 
should  not  present  any  intrinsic  hurdle. 

In  the  previous  section,  wc  have  introduced  the 
notion  of  the  principal  part  of  the  limit  of  the 
quantization  process  Eq.  (30).  In  the  initial  stages 
of  our  investigation  we  had  had  the  idea  that  by 
considering  the  constant  y  to  be  in  fact  a  slowly 
varying  function  of  x,  wc  could  find  the  actual  solu¬ 
tion  by  solving  for  7(1).  This  turned  out  to  be  a 
fruitful  pursuit  nnd  in  this  fashion  we  were  able 
to  find  excellent  numerical  replicas  of  the  solutions 
to  the.  differential  equation  without  any  tiouble 
except  in  the  neighborhood  of  a  turning  point.  This 
difficulty  had  to  do  with  choosing  a  trial  function 
for  the  logarithmic  derivative  of  the  particular  form 
(30);  the  problems  vanish  by  choosing  a  simpler 
starting  point. 

The  essential  ingredient  is  to  find  a  slowly  varying 
function  which  can  characterize  the  solution.  One 
possibility  that  we  have  just  discussed  is  to  use  a 
component  of  the  argument  of  the  tangent  function, 
that  is,  some  clement  dcsc'ibing  the  phase.  In  quan¬ 
tum  mechanics,  the  scattering  of  a  particle  by  a 
spherically  symm  trie  potential  can  be  described  in 
terms  of  phase  shifts.  By  introducing  a  nonlinear 
first-order  ordinary  differential  equation,  it  is  pos¬ 
sible  to  find  the  relevant  information  without  solving 
the  sccond-ordcr  Schrodinger  equation.  This  ap¬ 
proach,  introduced  by  Morse  and  Allis,'  and  more 


recently  discussed  by  Calogeio"  and  Ja*vy  and 
Keller,"  is  a  powerful  technique  for  finding  the  phase 
shift.  However,  these  authors  have  been  content 
with  specialized  information  and  have  overlooked 
the  possibilities  inherent  in  this  approach  which  can 
lead  to  solutions  of  the  sccond-ordcr  differential 
equation. 

Consider  the  function  <I>(x)  defined  by 

T(z)  =  -  k„  tun  (*,$(x)],  (32) 

where  kt  is  some  convenient,  but  arbitrary  reference 
value,  for  example,  the  mean  value  of  k*(x)  over 
some  interval.  Note  that  even  if  F(x)  behaves 
wildly,  the  corresponding  4>(x)  will  reflect  a  much 
calmer  character.  Since  F(x)  is  specified  once  4>(x) 
is  known,  we  can  ask  what  equation  must  4>(x)  obey 
if  r(x)  is  to  be  a  solution  of  the  Bicatti  equation. 
The  answer  is  found  by  straightforward  substitution 
into  (10),  and  after  simplification  we  find  that  4>(x) 
must  satisfy 

d*/dx  -  (A(x)  -  1]  cos’  (fr0«J>)  +  1,  (33) 

where  the  function  A'(x)  is 

A'(x)  *  k\x)/kl 

The  simple  formula  (33)  provides  all  the  information 
necessary  to  specify  a  precise  solution  of  the  second- 
order  differential  equation  -,.ith  virtually  no  com¬ 
plication  for  two  important  reasons.  First,  if  k2(x) 
is  bounded  then  the  derivative  i1'  is  always  bounded. 
Second,  this  formula  only  involves  K(x)  =  *’(x)/fcj 
so  that  imaginary  quantities  will  never  be  introduced 
if  tlie  original  differential  equation’s  coefficients  are 
real.  As  a  result,  it  is  a  trivial  matter  to  solve  equa¬ 
tion  (33)  by  an  iterative  process.  That  is,  given 
$,  =  #(x0)  at  x  «*  x«,  then  4>(x)  is  given  by  the 
expansion 

<t>(x)  =  4-  <f>iAr  +  i<J>r(Ax)’  ■+•••• 

+  ±  *r(Ax)'  +  •  •  ■  ,  (34) 

n  1 

where  the  higher-order  derivatives  arc  easily  cal¬ 
culated  from  (33), 

$"  =  -  1)  sin  (2M-)  +  A"  cos’  ( k ,$) 

-  —  M>"(A'  -  1)  sin  (2A0$)  -  2 *.*' 

X  {A"  sin  (2fto4>)  +  UK  -  1) 

X  cos  (2Ayt>)]  +  K"  cos’  (*<,$) 

•  * 

$'•’  -  etc.  (35) 

’*  F.  Calogero,  Nuovo  Cinirnto,  27.  261  (1063). 

»  B.  R.  Levy  nnd  i.  B.  Keller,  J.  Msth  Phy.  4,  54  (1953). 
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In  other  words  we  have  a  recurrence  relation 
expressing  the  higher-order  derivatives  in  terms  of 
known  lower-order  ones.  Once  we  are  given  <J'0|  it  is 
no  problem  to  determine  what  4>,  —  *f*(i ,)  is  at 
some  nearby  point  x,  to  any  degree  of  accuracy 
,by  truncating  the  expansion  (34).  A  repetition  of  the 
process  allows  one  to  calculate  <h:  *=  *h(jf)  and  so  on. 
No  doubt  some  purists  will  argue  that  this  is  not 
a  closed  form  representation  of  the  solution.  While 
this  is  true,  it  is  in  large  part  a  meaningless  objec¬ 
tion.  Any  number  of  closed  form  expressions,  c,g. 
contour  integrals,  are  very  elegant  and  eoncisc  but 
quite  intractable  when  they  are  interrogated  for 
numerical  data.  The  important  question  to  be 
answered  is  how  readily  can  desired  information  be 
obtained.  Before  we  can  give  any  demonstrations 
of  the  utility  of  Eq.  (33),  we  must  digress  and 
examine  the  intractability  of  a  closed  form  solution, 
namely  Eq.  (3). 

We  have  implied  that  once  the  logarithmic  deriva¬ 
tive  T(x)  is  known,  we  could  reconstruct  the  solution 
ti(x)  from  the  integral  (3).  While  theoretically  this 
formula  should  provide  the  desired  answers,  the 
practical  situation  is  otherwise  since  the  integral 
diverges  at  the  poles  of  r(x).  If  u'(x)  is  bounded, 
then  any  such  infinity  in  T(x)  is  just  another  way 
of  saying  that  u(x)  has  a  zero  at  that  point,  but 
nonetheless  numerical  integration  of  Eq.  (3)  to  ob¬ 
tain  u(x)  is  to  a  great  extent  out  of  the  ques  ion. 
In  other  words,  even  with  an  explicit  T(x)  at  our 
disposal,  we  may  have  all  sorts  of  grief  in  translating 
this  knowledge  into  numbers  describing  «(x)  by  use 
of  (3)— so  much  for  closed  form  expressions. 

The  scheme  we  prefer  to  adopt  neglects  formula 
(3)  and  instead  interpolates  u(x)  by  a  polynomial 
at  each  stage  of  the  iterative  process.  At  some  x,, 
suppose  we  know  I\  and  an  initial  value  at 
x..,  >  x,  we  arc  given  T..,.  These  three  pieces 
of  information  determine  a  u:iiquc  sceond-ordcr 
polynomial  which  must  assume  the  value 


tt.-i 


u. 


2  +  r,(xiM  -  x,) 

2  ■“  F,*,(x#*i  Xi) 


(36) 


at  x  =  x.*,  (the  derivation  of  this  formula  is  given 
in  the  Appendix).  The  ability  of  this  rational  func¬ 
tion  to  predict  the  future  value  u.>,  in  terms  r(, 
and  r,.|  is  rather  remarkable.  For  example,  with 
an  increment  r,»,  -  x(  =  0.1,  Eq,  (3G)  reproduces 
sin  x  from  its  logarithmic  derivative  cot  x  to  within 
four  significant  figures.  In  addition,  Eq.  (3G)  auto¬ 
matically  yields  a  zero  of  u(x)  should  T(x)  have  a 
pole  at  a  point  where  u(x)  is  regular,  so  that  the 
role  of  the  singularities  of  T  is  correctly  reproduced. 


With  these  preliminaries,  we  can  find  numerical 
solutions  of  seen, id-order  differential  equations  to 
nny  desired  degree  of  accuracy.  The  iterative  scheme 
which  carries  the  solution  from  one  point  to  another 
is  to  truncate  the  Taylor  expansion  (34)  to  a  poly¬ 
nomial  of  suitable  order  for  the  increment  chosen. 
A  numerical  example  may  serve  to  illustrate  the 
procedure.  The  differential  equation 

d'u/dx1  *f  (a*  —  C/x’)  u  ■»  0  (37 ) 

is  distinguished  by  having  a  particularly  simple  pair 
of  solutions 

w,  =>  (3/ox)  cos  ax  +  (1  —  (3/a*x’)J  sin  ox 


+  (l  ~  »>n  (a*  +  jj)*  (38) 

Since  (37)  has  a  turning  piont  at  x  =  6 V’a  as  well 
as  a  singularity  at  the  origin  it  poses  a  rather 
severe  computational  challenge.  Wc  have  used  this 
equation  as  a  '  st  of  our  procedure  for  various 
choices  of  the  parameter  a.  The  larger  this  parameter, 
the  more  rapid  are  the  oscillations  of  the  solution 
for  jxj  >  GVa,  and  the  more  likely  arc  errors  to 
grow  in  any  computational  scheme.  By  use  of  a 
quadratic  truncation  of  (35)  wc  have  found  that 
the  formulas  (35)  and  (36)  reproduce  the  solution 
very  satisfactorily.  As  mentioned,  the  errors  get 
worse  with  increasing  a  and  Fig.  2,  illustrates  one 
of  the  poorer  cases  wc  have  considered  for  which 
the  parameter  a  had  the  value  3.5  and  G*/a  ■=  0.70. 
The  increment  chosen  was  Ax  =  0.0142,  a  rather 
odd  number  which  resulted  from  choosing  an  in¬ 
crement  inversely  proportional  to  a  and  equal  to 
0.1  if  a  =  0.5.  At  x  *=  —10,  the  initial  logarithmic 
derivative  was  adjusted  to  match  u2,  the  solution 
with  a  singularity  at  the  origin.  From  this  initial 
value  a  solution  was  computed  in  the  direction  of 
increasing  x  using  a  quadratic  formula  for  4>(x)  to 
compute  T(x).  This  was  followed  by  use  of  Eq.  (30) 
to  determine  the  solution.  The  calculation  and  the 
actual  solution  agree  so  well  that  it  is  not  worth¬ 
while  to  draw  curves  of  each.  Instead,  the  absolute 
error  curve  is  drawn  with  a  magnified  scale.  As 
might  bo  expected  in  any  initial-value  calculation, 
the  errors  tend  to  grow  as  the  computation  pro¬ 
gresses.  However,  except  in  the  neighborhood  of  the 
singularity  we  note  that  the  amplitude  of  the  error 
increases  but  linearly,  and  furthermore  it  oscillates 
about  zero.  The  divergence  at  the  origin  is  more 
apparent  than  real  since  the  percentage  error  remains 
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PER  linear  differential  equations 


boun.loa  By  choosing  a  smaller  increment  or  addins 
L  *al  (orln  to  the  expansion  for  'l’(x),  the 

r  can  be  reduced  by  an  additional  order  of 
magnitude.  Wc  might  add  that  the  eomr^^ 

.re  sufficiently  simple  and  proccrn  so  fast  that  they 
can  §bc  performed  with  but  the  aid  of  a  d^L  enl- 
cuLor  Kven  on  a  relatively  slow  IBM  1  HO  com¬ 
puter  the  calculations  were  performed  as  fast 
the  machine  could  print  out  the  answers. 

Thp  preceding  example  illustrates  the  case  with 
wUch  a sccond-onler  equation  can  be  solved  wuhout 
any  complications  or  need  ot  any  sophistics, 
numerical  techniques.  Of  course,  certain  tacks  of 
thTtrade  can  be  very  beneficial  >n  reducing  errors. 

For  example  we  recommend  use  of  a  variable  n 
crement  which  is  smallr  l  the  initial  atages  of  1  c 
S Union  to  reduce  a  buildup  of  error,  In  add. mu 
the  increment  ....  -  ..  ^ould  be  scaled  so  that 
*  _  *4  docs  nut  exceed  some  preassigned l bound 

In’ this  fashion,  a  computer  program ,  c  an  be  both 
sDcedcd  up  and  made  more  reliable  by  usirg 
Eg  increment  whose  size  vanes  inversely  with 

the  size  of  the  change  in  *(*)  from  x,  to  *i... 
pLlly  we  note  that  if  the  logarithmic  derivative 
l  o(  Eq  (6)  is  known,  then  ‘-he  logarithmic  der.va-  ( 
S;e  r^0f  aquation  (1)  »hlcfc  i«pW  U»  , 

is  given  by 

r„  -  r.  -  MW  ^ 

so  that  without  any  further  ado,  the  solution  of 

^  a  lwt  remark  te  should  like  to  point  out  that 
the  results  of  this  section  complement 
supersede  the  discrete  analysis  of  Sec.  II  for  the 
Z  of  a  differential  equation  with  analytic  coeffi¬ 
cients.  For  whereas  this  section  is  concerned  with 
accurate  representations  which  need  a  relatively  fine 
interval  for  their  calculation,  the  analysis  m  Sec.  II 
yields  good  approximation  for  coarser  interxa  . 
While  this  may  be  of  no  matter  when  an  equation 
U  to  be  solved  but  once,  the  distinction  bccom 
‘  important  when  the  differential  equation  must 
wived  many  time,  to,  w.ue  wl  o(  ixrtoha  »n» 
in  its  coefficients.  The  formulas  in  Sec.  II  furnish 
a  quick  method  of  obtaining  insight  into  the  effect 
:fr, ^number  of  variations  and  locating  regions 
demanding  further  investigation.  Following  this  lx- 
nloratioi,  the  formulas  of  See.  Ill  can  chart  these 

unknown  domains  To  cite  a  nautical  analogy :  it  is 
the  difference  between  using  dead  reckoning  and  an 
inertial  guidance  system.  One  is  cheap  »n  0 
adequate,  the  other  is  more  expensive  but  arbitrarily 

accurate. 
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Flo  2.  The  upper  curve  is  the  eolation  Ui(i)  o(  Fq  bJJ). 
and  the  lower  curve  is  the  computational  error  of  the  numerical 
procedure  described  in  the  text.  Note  that  there  is  a  1U.I 
ratio  between  the  two  scales. 
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appendix 

In  this  Appendix  wc  shall  derive  the  formula  (36) 
in  Sec.  III.  From  the  equation  (9)  we  have 

uT  =  u'.  (4°) 

|  With  no  loss  of  generality  suppose  x.  =  0.  Let  us 
.  assume  that  the  function  u(x)  can  be  approximated 
3  by  some  unknown  quadratic  expression  which  can 
.  be  chosen  in  the  completely  general  form 

°  u(x)  =  *x(x  -  A)  +  0(x  -  A)(x  +  A) 
a  +  yx(x  +  A),  (41) 

y  where,  o,  0,  7  are  constants  to  be  determined  and 
A  is  the  length  of  the  increment. 
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Substituting  (41)  into  (40)  wo  find  that 
lx' (a  +  0  +  y)  +  *A(-a  +  y)  -  0A*]r 


- 

2x(a  +  A  +  y)  +  A(  —  a  -f-  y). 

(42) 

Let  u,, 

ii.  aud  r»,  r#,  r_  denote  the  values 

of  u  nud 

r 

nt  the  points  x  ■=■  A,  x  “  0,  * 

-  -A, 

respectively 

.  Then,  from  (42)  we  obtain 

u.T* 

M 

(a  +  20  +  3y)A  -  2yA*r*, 

u,r. 

MM 

(-«  +  y)A  -  -0A‘r* 

(43) 

u.r. 

- 

-(3a  +  20  +  y)A  -  2A*ar.. 

The  last  two  equations  imply  that 

w-  _  J2a  ^  2  -  Af, 
u.  "  “4  "24-Ar. 

or 


(44) 


u,  -  «-(2  +  AF_)/(2  -  AT.)  (45) 

By  similar  reasoning  more  lophisticatcd  formulas 
can  be  derived  whose  accuracy  is  of  arbitrarily 
high  order. 
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Travel  Time  and  Phase  Shift’ 
J.  Kane1 

California  Imtilulc  of  Technolog  j ,  Pawfena 


Abstract.  U  any  harmonic  wave  propagate:;  over  a  path  with  a  nommiform  velocity,  then, 
owing  to  second -order  effects,  pha«c  shift  data  will  not  lx>  proportional  to  travel  time.  The 
travel  time  will  alway  rccni  to  he  larger  than  it  actually  is  owing  to  increased  phase  shift 
resulting  from  diffraction  effects.  This  dilation  can  he  described  in  a  rather  elementary  manner 
in  temia  of  a  nonlinear  first-order  differential  equation.  An  integration  of  this  equation  gives 
an  explicit  formula  for  the  I'  screpamy  between  actual  travel  time  and  the  travel  lime  com¬ 
puted  i.'om  phase  shift  data. 


Analysts.  An  often  drawn  eonelus  n  from 
Fourier  analysis  is  the  interchangeability  of 
phase  information  in  the  harmonic  domain  and 
travel  time  in  the  real  domain.  While  this  con¬ 
nection  can  be  firmly  established  for  uniform 
wave  propagation,  its  applicability  to  propaga¬ 
tion  in  a  perti  rbed  continuum  is  'css  sccut  In 
this  no‘e  »e  will  consider  this  question  from  a 
point  of  view  which  .".Slows  us  to  derive  a 
,'innula  for  the  phase  shifts  of  harmonic  waves 
traversing  an  inhomogeneous  vt.^city  profile. 
This  permits  us  to  make  an  explicit  comparison 
between  travel  time  and  phase  shift,  its  pre¬ 
sumed  harmonic  image.  The  procedure  we  fol¬ 
low  Is  to  introduce  a  nonlinear  first-order  equa¬ 
tion  Tor  the  ph.ae  variation.  In  terms  of  this 
eauation  we  can  compute  the  phase  shift  be¬ 
tween  two  (mints  as  an  integral  involving  the 
velocity  profile  without  the  need  of  solving  the 
wave  equation  This  provides  an  explicit  basis 
for  a  comparison  between  phase  shift  data  and 
travel  time. 

For  reference,  we  first  cite  the  formula  giving 
the  travel  time  Tm  of  a  mass  particle  moving 
with  variable  velocity  r(x)  between  the  two 
points  x,  and  i,,  n  .ncly 


From  the  point  of  view  of  causality,  wt  ex¬ 
pect  that  this  formula  would  also  express  the 


*  Cmitribu'ion  1303,  Diiirion  of  Geological  Sci¬ 
ences,  Cahtorria  Institute  of  Technology,  Pasa¬ 
dena. 

*  On  academic  •cove  from  the  University  of 
Rhode  Island. 


travel  time  of  a  harmonic  wave  propagating  in 
s.  pertnrlicd  continuum.  The  question  wc  ask  is 
whether  (1)  is  consistent  with  phase  shift  data. 
This  query  has  significance,  tor  wave  motion 
has  features  which  arc  foreign  to  the  motion  of 
a  mass  point — for  example,  continuous  reflection 
and  transmission.  For  a  harmonic  wave,  x 
pcct  that  the  cumulative  effect  of  diffractio. 
phenomena  will  be  included  in  its  phase  varia¬ 
tion.  In  short,  for  nonuniform  wave  propaga'ion 
wc  anticipate  that  the  phase  shift  between  two 
points  will  have  scattering  effects  corrupting  the 
travel-time  information.  Wc  want  to  find  a 
formula  that  separates  these  two  effects  so 
that  phase  information  can  be  corrected  to  yield 
the  travel  time. 

Uniform  steady-state  wave  motion  is  described 
by  tlie  wave  equation 

d’u/dx1  +  A'o’t'M  —  0  (2) 

where  the  constant  A-«  is  inversely  proportional 
to  the  velocity  of  propagation 

k0  =  ai/f0  (3) 

From  a  knowledge  of  the  piano  wave  solution 
it  is  possible  to  measure  the  travel  time 
T.  between  x,  and  x.  of  a  harmonic  wave  by 
dividing  the  phase  diffc  ence  A<1>  by  the  fre¬ 
quency  <i» : 

T .  --  '\i  =  (x2  —  X,)  Co  (4) 

As  long  as  the*  velocity  is  uniform,  the  solutions 
of  (2)  are  obtain-  J  by  inspection,  and  :his  tech¬ 
nique  offers  no  difficulty.  However,  complexities 
arise  if  the  wave  traverses  an  inbomogenous 
path  "long  which  the  velocity  is  no  longer  uni- 
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form,  for  then  the  wave  equation  w  not  amen¬ 
able  to  direct  solution. 

Let  u«  write  the  variable  velocity  v(z)  ns 

r(x)  -  .•„/! I  +  ‘Ml  <0 

where 

“  [v#/r(x)]  -  (6) 

describes  the  local  velocity  perturbation.  To 
compute  the  travel  time  T,.  we  need  to  calcu¬ 
late  the  phase  shift  of  the  traveling  wave  solu¬ 
tion  of  the  pcrturlxnl  or.  tation 

d'u'djc'  4*  £o'(  1  +  «(x)l’u(x)  =  0  (7) 

between  the  points  /,  and  x,  as  compared  with 
the  phase  shift  A<1>  =  A.(x,— x,)  for  the  uniform 
case.  For  nonnnifonn  wave  propagation  it  is  not 
possible  to  identify  a  traveling  wave  without 
ambiguity,  for  there  can  be  no  unique  distine 
tion  between  the  forward  and  scattered  waves 
[SchclLunoff,  1951].  Hence,  the  only  specifica¬ 
tion  coneem  ng  the  behavior  of  u(x,  t)  that  we 
»hall  make  is  that  it  becomes  Re  e"**"*"  as 
«  -»  0.  For  clarity,  we  shall  concentrate  on  real 
wave  motions;  that  is,  we  discuss  fields  in  which 

u(x,  t)  ~  cos  (fc,,x  —  «*><)  (8) 

Note  that  this  forces  us  to  include  the  effects  of 
timevar.  ioninufx,  I). 

To  obtain  the  phase  we  first  introduce  the 
logarithmic  derivative 

r(x,  t)  *=  (du/dr)/u  (9) 

and  in  terms  of  this  quantity,  (7)  becomes  the 
Ricatti  equation 

ar/dr  +  r  +  *,’d  +  «)’  =  o  (io) 

If  «(i)  vanishes,  the  solution  of  (10)  subject  to 
(8)  is 

r'x,  0  “  —  *0tan  {V  —  «<]  (II) 
If  «(r)  is  nonzero,  we  seek  a  solution  in  the  per¬ 
turbed  form: 

r(a,  I)  «  -fc>  tan  jfc0(x  +  ?(*)]  -  w(j  (12) 

wlic.’  h\y(x)  descrilrcs  the  local  phase  correc¬ 
tion  owing  to  a  nonnniform  velocity.  That  is, 
the  total  phase  shift  A<f*  will  be 

A#  -  -  *(xt)  4-  xi  —  *i)  (13) 


In  terms  of  y(z)  as  an  unknown,  the  Ricatti 
equation  becomes 

d? /di  **  (2 1  +  »’) 

•cos1  |<r0(x  +  *>(/))  —  ut\  (14) 

after  substitution  of  the  trial  form  (12)  and 
some  arithmetic.  If  we  integrate  this  expression, 
we  obtain  an  intrgril  equation, 

Alt)  ~  A* i)  *=  f  (2»(t)  +  «’({)) 

•'ll 

•cos’  |t0[i  +  vr(i) J  —  <*d|  rf(  (15) 

for  the  dctrrinin.it ion  of  ^(x),  sinre  it  appears 
under  the  integral  sign.  Note  that  the  net  cor¬ 
rection  to  the  phase  shift  f  (x.)  —  j>(x,)  is  time- 
dependent.  This  situation  is  forced  bv  our  in¬ 
ability  to  separate  the  forward  wave  from  the 
multiply  reflected  wavelets  it  excite*.  Only  i. 
uniform  wave  motion  is  the  phase  time-inde¬ 
pendent,  for  then  the  real  and  imaginary  parts 
of  the  exponential  are  identical  apart 

from  a  translation.  In  this  special  ease,  the 
shape  of  the  wave  docs  not  change  ns  we  pass 
through  a  period.  For  nonuniform  trove  propa¬ 
gation  the  shape  of  the  wave  changes  during  one 
cycle,  and  for  this  reason,  we  cannot  speak  in 
general  of  a  distinct  phase  shift  between  two 
points  owing  to  its  time-variation.  However,  ex¬ 
perimentalists  often  measure  phase  shift  aver¬ 
aged  over  many  cycles.  If  this  is  done,  the  cos’ 
{}  term  in  the  integral  equation  (13)  behaves 
like  a  numerical  weight  factor  of  ],  and  we  ob¬ 
tain  the  explicit  relation 

(A^)  ~  Ait))  =  f  («(£)  +  $«’(£)]  >f£  (16) 

where  the  angle- brackets  0  indicat"  an  average 
phase  difference.  Alternatively,  (16)  describes 
what  would  be  measured  if  x,— x,  is  large  in 
comparison  with  a  wavelength.  Thus,  the  travel 
time  r,  for  harmonic  waves  in  terms  of  the 
average  phase  shift  is 

Tw  =  (x,  -  x,)  +  (Ah)  -  sp(x,)>  (17) 

or 

r.  -  f‘  { 1  +  <({)  +  *<’({)}  d(  (is) 

as  compared  with  (1).  This  expression  can  be 
rewritten  as 
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t.  -  r.  +  .4  '>9) 

where  the  anomaly  /4  lx-twecn  T .  mui  . is 
wn  ro  be 


and  we  always  lux  e 

T.  >  Tm  (21) 

and  equality  holds  if  ami  only  if  i-(x)  ra  iy  We 
interpret  tins  anomaly  .4  as  a  result  of  the  in¬ 
creased  phase  shift  due  to  diffraction. 

Example.  The  implication  of  the  preceding 
remarks  is  ix-rhans  best  illustrated  by  an  ex¬ 
ample.  Consider  a  surface  wave  propagating 
along  a  path  whose  layer  thickness  varies  be¬ 
tween  A...  and  A,,,,  say.  For  a  uniform  layer 
we  can  find  a  period  equation  which  rrl  ites  the 
phase  velocity  of  any  harnio:  ;  surface  wave 
and  the  height  A.  The  question  we  ask  is  what 
will  the  average  height  (A)  scent  to  be  if  we  cor¬ 
relate  phase  shift  data  over  large  distances  with 
phase  velocity. 

In  its  rigorous  formulation  this  poses  an  in¬ 
tractable  boundary  value  problem,  lint,  if  the 
variation  of  A  is  small  over  distances  of  the 
order  of  a  wavelength,  then,  the  propagation  of 
the  disturbance  is  described  by  the  one-dimen¬ 
sional  equation  (7)  where  i- (x)  is  the  local 
phase  velocity.  Zither  the  local  phase  velocity 
t-(x)  or  the  height  A(x)  determines  the  other 
owing  to  the  |ieriod  equation.  ror  simplicity,  let 
us  assume  that  their  interrelationship  is  a  linear 
one,  although  the  argument  we  shall  present 
does  not  de|>end  upon  this  fact.  If  linearity 
holds,  the  arithmetic  mean  of  the  velocity  t>(x) 
w-ould  be  proportional  the  average  thickness  (A). 
How  can  we  extract  the  arithmetic  mean  from 
measured  phase  data?  We  write 

c(x)  -  r9[i  -f  i(x)j  (22) 

where  r.  is  the  arithmetic  mean  of  t-(x)  and 


4{x)  is  related  to  <(/)  by  the  relation 

«(/)  =  -4(x)'(!  +  4(/)|  (2d) 

Since  r,  is  the  arithmetic  ttv  in,  we  have 

r  m  o  (2i) 

but  the  mean-square  value 

<A\'x)>  =  - —  f’  «’({)  df  (?h) 

x»  x,  J,, 

is  nonzero,  For  a  fixed  distance,  velocity  is  in¬ 
versely  proportional  to  travel  time.  From  (1) 
we  have  the  second-order  equation 

(«')-  “  <•->/[!  +  <i’(-r))l  (26) 

as  the  effective  velocity  of  a  mass  point  in 
traversing  the  distance  from  r,  to  x,  If  we  in¬ 
terpret  (t’>.  as  the  velocity  we  measure,  (1  + 
<5*(x) >)  is  the  eorrertion  faetor  we  need  to  re¬ 
normalize  it  to  the  arithmetic  mean.  However, 
the  inequality  (21)  indict  '£  that  this  correction 
factor  is  too  small  for  a  wave  propagating 
through  an  inhomogeneous  velocity  profil-  owing 
to  scattering  effects.  In  'act,  for  a  wave,  we 
easily  find  that  the  effect  .vc  velocity  is  related 
to  i’,  by  the  second -order  relation 

<V>.  =  !u/(l  +  J(4’(x))]  (27) 

In  other  words,  the  effect  of  scattering  is  to  in¬ 
crease  the  presumed  correction  factor  by  507c. 
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Fiber  Optics  and  Strain  Interferometry 
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Abstract — By  use  of  glsss  fiber*  as  d  'ectric  waveguides,  lnter- 
feromeiric  measurements  of  light  along  t.-  .  propagation  path  can  be 
used  as  a  new  type  of  strain  gauge.  Two  novel  strain  interferometer* 
are  proposed:  one  which  uses  the  phenomenon  of  crosstalk  between 
adjacent  light  pipes,  and  another  which  Is  a  waveguide  version  of  the 
Michelson  interferometer.  Both  designs  are  distinguished  by  the 
possibility  of  digital  readout  and  a  basic,  simplicity.  The  theory  of 
the  instruments  is  described,  and  explicit  estimates  for  the  Strain 
sensitivity  are  calculated.  The  numerical  data  show  that,  with 
digital  readout,  strains  of  the  order  of  10  '  can  be  resolved  with  a 
basic  instrument  ten  meteia  long,  and  strains  of  the  order  of  10~* 
can  be  resolved  by  a  hundred-meter  device.  With  the  aid  of  various 
refinements,  the  sensitivity  can  be  enhanced  by  at  least  another 
order  of  magnitude. 


1.  Introduction 

GROUND  MOTION  can  be  measured  with  rea¬ 
sonable  accuracy  by  a  wide  variety  of  seismom¬ 
eters.  However,  other  than  the  BeniofT  strain 
seismometer  or  an  installation  which  uses  multiple  beam 
interferometry  between  piers  in  a  seismic  vault,  there 
are  few  instruments  that  arc  capable  of  recording  the 
slow  accumulation  of  tectonic  strain.  The  demands  of 
the  projected  earthquake  prediction  program  indicate  a 
need  for  a  network  of  such  strain  gauges  which  combine 
modest  cost  with  adequate  performance.  In  this  study, 
wc  investigate  the  possibility  that  an  amalgam  of  wave¬ 
guide  and  optical  techniques  may  be  useful  in  realizing 
these  goals. 

Seismologists  have  long  used  surface-wave  dispersion 
data  to  gauge  the  layering  of  crustal  structures.  How¬ 
ever,  this  nonlinear  interplay  of  dispersion  and  depth 
is  a  feature  of  wave  motion,  in  general,  and  not  con¬ 
fined  to  the  theory  of  elastic  waves.  Indeed,  optical 
surface  waves  have  been  observed'  !  and  their  features 
studied.1'4  Owing  to  the  minute  character  of  optical  sur¬ 
face  waves,  interferometric  procedures  based  upon  their 
interaction  with  matter  can  supply  rather  accurate  in¬ 
formation  about  small  displacements  or  strains.  In  this 
analysis  we  propose  an  application  of  optical  surface 
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waves  propagating  along  dielectric  fibers  as  a  new  type 
of  displacement  or  strain  gauge  with  digital  readout. 

The  basic  instrument  can  be  thought  of  as  a  long 
thin  prism  with,  for  example,  a  square  cross  section  in 
the  quiescent  state.  Under  the  at  .ion  of  forces  normal 
to  the  axis  of  the  prism,  tite  square  cross  section  w  ill  be 
deformed.  Depending  upon  the  stress  distribution,  the 
strained  cross  section  can  assume  a  variety  of  forms: 
compressed  square,  rhombus,  etc.  The  distorted  geom¬ 
etry  can  he  characterized  by  measuring  the  perturbed 
distances  between  the  edges  of  the  prism  and  comparing 
them  with  those  of  the  equilibrium  state.  At  optical  fre¬ 
quencies,  waveguides  can  be  very  small,  and  we  can 
imagine  that  the  four  edges  arc,  in  b  et,  miniature  wave¬ 
guides.  If  the  waveguides  are  coupled,  then,  they  can 
interact  in  such  a  fashion  as  to  respond  to  the  varying 
distances  betw  een  them.  Wc  can  imagine  that  these  four 
waveguides  arc  relatively  dense  dielectric  fibers  cladded 
in  a  transparent  matrix  to  produce  a  multicomponent 
glasi  hie.  i  nis  configuration  would  then  be  imbedded 
in  some  elastic  medium  to  couple  it  to  that  structure. 
For  example,  it  could  be  buried  in  the  earth,  or  im¬ 
bedded  in  a  concrete  dam  or  foundation.  Once  inserted, 
the  average  distortion  of  the  transverse-  cross  section  can 
be  monitored  by  studying  the  variation  of  the  propaga¬ 
tion  characteristics  along  these  coupled  microscopic 
waveguides.  The  purpose  of  this  paper  is  to  describe 
how  the  interaction  of  coupled  waveguides  can  be  used 
as  a  metrological  tool  for  the  determination  of  sinc.ll 
displacements  or  strains. 

The  principle  wc  propose  using  to  achieve  our  goals 
is  an  application  of  the  phenomenon  of  crosstalk  which 
has  been  observed  in  the  optical  range  of  couplcd-modc 
propagation.  This  is  an  exchange  of  energy  between 
adjacent  light  pipes  which  lias  been  analyzed  by  Snitzcr* 
who  had  in  mind  the  problem  of  controlling  moding  in 
fiber  optic  lasers.4  In  this  regard,  crosstalk  is  a  problem 
and  the  motivation  is  to  reduce  this  nuisance.  Our  in¬ 
terest,  however,  is  to  propose  its  exploitation.  For  this 
purpose  wc  review  the  theory  of  this  effect,  obtain  esti¬ 
mates  for  the  sensitivity  that  can  be  obtained,  and 
propose  explicit  designs  for  working  iastruments. 

'The  basic  building  blocks  of  the  device  are  sketched 
in  Fig.  1.  Two  dielectric  w  aveguides  arc  separated  from 
one  another  by  a  distance  d.  Each  waveguide  is  of  di¬ 
mension  h  and  is  a  dense  dielectric  of  index  imbedded 
between  rarer  media  « i  and  tit.  The  geometry  of  the 
device  can  consist  of  either  a  pair  of  cylindrical  fibers,  a 
pair  of  plane  parallel  layers,  or  any  other  pair  of  parallel 
waveguides,  provided  that  the  structure  ha*  a  plane  of 
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symmetry  l»i*a cting  the  distance  d  separating  the  com¬ 
ponent  waveguide*.  Since  the  operation  of  the  device 
depends  primarily  upon  symmetry  considerations  and 
not  especially  it  iron  the  character  of  the  component 
waveguides,  we  shall  provide  explicit  calculations  only 
for  plane  parallel  structures  for  ease  of  analysis.  How¬ 
ever.  we  stress  that  the  theory  is  easily  generalized  to 
describe  other  structures.  An  mule-standing  of  the  de¬ 
vice's  function  depends  upon  a  knowledge  of  guided 
wave  propagation  in  this  multilayer  system.  That  is,  wc 
sock  the  dependence  of  the  system’s  behavior  upon  the 
parameter  d  which  is  the  quantity  wc  should  like  to 
measure. 


Fig.  t.  The  geometry  of  a  symmetrical  two-fiber  system. 

II.  Theory 

To  analyze  this  multilayer  system  we  need  to  find 
those  solutions  of  Maxwell’s  equations  which  satisfy  the 
boundary  conditions,  and  which  represent  propagating 
modes  for  our  choice  of  geometry.  If  we  denote  the 
spatial  part  of  the  electric  vector  as  E  and  the  spatial 
part  of  the  magnetic  vector  as  H,  then  in  homogeneous, 
isotropic  space,  Maxwell’s  equations  assume  th  form 

curl  E  =  ikJH  (1) 

curl  H  «  —  iktHj-E  (2) 

where  a  time  factor  c-'"'  has  been  suppressed.  The  sym¬ 
bol  k9  is  the  propagation  constant  of  free  space 

2»/A»  —  u/c]  (3) 

and 

kj  *=  kt»j  (4) 

where  we  identify  X*  and  c  as  the  wavelength  and  speed 
of  light  in  vacuum,  and  n,  is  the  dielectric  constant  in 
thejth  mHium.  For  convenience,  the  magnetic  permea¬ 
bility  n  is  assumed  to  be  unity  for  all  media,  and  we 
neglect  loss  terms  proportional  to  conductivity. 

If  the  wave  disturbance  travels  along  the  x-axis  in  an 
isotropic  medium  without  dependence  on  the  Cartesian 
coordinate  s,  then  all  derivatives  with  respect  to  z  van¬ 
ish.  One  can  show  from  (1)  and  (2)  that  only  the  two 
transverse  components  U(x,  y)  =  E,  and  T(x,  y)  «  H,  are 
required  to  specify  all  six  components  of  the  electromag¬ 
netic  field.  These  six  components  separate  into  two  de¬ 
coupled  groups  of  three  each.  Thus, 


+  i  au  -1  at 

Te  '  ™  V~  ’  //»  = - —— 

lL-  1  dy  ik  ^  ds  (5) 

fit  '■  fit  *=  fi,  0, 

determine  Et  and  E „  f.ou.  U(x,  y)  for  TE  propagation, 

and 

T***  *  tt— -  -  -7- - 

«*»«/  *y  i*#V  dx  (6) 

F.t  =  //,  -  11^  »  0, 

specify  H,  and  Hs  in  terms  of  K(x,  y)  for  7" M excitation. 
Both  V  and  V  satisfy  the  same  reduced  wave  equation, 
namely; 

(V*  4  kflV,  V  -  0,  (7) 

in  the 7th  medium  which  has  index  of  refraction  Hj.  The 
field  within  the  7th  medium  must  satisfy  this  wave 
equation  and  the  requirement  that  the  tangential  com¬ 
ponents  of  E  and  H  be  continuous  across  any  interface 
separating  two  dielectrics.  For  TE  modes,  this  means 
that  U  and  dU/dy  must  be  continuous  across  any 
boundary;  fo.  TM  modes,  V  and  (1/kfldV/dy  must  be 
continuous. 

The  nature  oi  the  problem  to  be  solved  is  such  that 
only  discrete  types  of  propagation  are  permissible  along 
the  system.  For  plane  layers,  the  electromagnetic  theory 
is  very  similar  to  that  of  Love  wave  propagation.  Since 
the  TEand  TM  problems  arc  identical  in  character,  wc 
shall  emphasize  the  analysis  for  the  TE  case  and  just 
quote  the  corresponding  results  for  the  TM  case, 
hor  a  traveling  mode  we  must  have 

U(x,  y)  «  «.v{y)e(*»v%  (8) 

and  then  the  problem  becomes  one  of  determining  those 
values  of  the  dimensionless  propagation  constant  N  for 
which  functions  «,v(y)  can  be  found  which  a/e  consistent 
w  ith  Maxwell’s  equations  and  the  boundary  conditions. 
At  each  interface,  wc  have  two  transition  conditions  to 
,bc  satisfied.  If  we  take  the  ratio  of  these  equations,  wc 
obtain  an  equation  whose  physical  significance  is  inde¬ 
pendent  of  considerations  of  the  field  amplitude.  For 
this  reason  we  introduce  the  logarithmic  derivative 

,  du  / 

rw  *  jJ  *  <»> 

which  is  essentially  an  impcdancc-typc  function.  For 
the  TE  case,  it  is  proportional  to  H,?E„  which  is  the 
normal  admittance  to  any  guide  wall  where  it  must  have 
the  same  value  011  either  side  of  the  boundary.  For  the 
TM  case,  the  ratio  17 ny  is  an  impedance  function  and 
plays  an  analogous  role.  In  the 7th  region  u.v(y)  must  be 
of  the  form 
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Mi  cos  (*.(«,’  -  A’*),,tjf] 

1  +H,  sin  [*.(«,*  -  iV  <  Hf 

"‘v(y)"  i,l#«p[+*.(A-*-« Z)1^]  ’  ' 

[  +  fl/cxp  [— i»(.V*  -  A"  >  w, 

and  we  must  solve  ait  algebraic  problem  to  determine 
the  correct  values  of  A,  and  Hf  so  that  the  transition 
conditions  across  each  interface  arc  matched  and  radia¬ 
tion  conditions  at  infinity  arc  satisfied.  Instead  of  pro¬ 
ceeding  directly,  computational  effort  can  decreased, 
and  insight  enhanced,  by  considering  the  operation  of 
the  system  as  a  coupling  of  the  J Turactcristics  of  the 
individual  waveguides.  We  analyze  the  behavior  of  an 
individual  layer  first,  and  then  couple  two  guides  to¬ 
gether  to  understand  the  response  of  the  composite  sys¬ 
tem. 

III.  Modai.  Analysis 

A.  The  Component  Guide:  Ah  isolated  Layer 

In  order  to  discuss  propagation  along  the  structure 
sketched  in  Fig.  1,  we  assume  that  the  /-variation  be¬ 
haves  like  a  standing  wave  within  the  layer,  and  an 
exponentially  damped  field  outside.  These  fields  will 
match  across  the  interfaces  in  such  a  way  that  the 
transition  conditions  are  satisfied  if  the  allowable  values 
of  the  propagation  constant  N  satisfy  the  characteristic 
equation.* 

tan  [(«,*  -  X*)utkji] 

(Ht*  -  N')'”(r,m  +  r»)  TE 

(m*  -  .v*)  -  r»-  r, 

for  the  TE  case  and 
tan  [(«,*  -  y*yi*kth] 

(«»’  -  .v*)1'i[(n,/n,)*r,«  +  («!./«, )*r,]  _  ,.„x 

■  ■  1  "  . — —  Tin  (12} 

(a,*  -  ,V»)  - 

for  the  TM  ca»e.  In  these  expressions,  the  indexed  sym¬ 
bols  are  given  explicitly  as 

r,-  -  VW’  -  «i* 


however,  the  field  outside  is  perturbed  in  such  a  fashion 
that  the  wall  admittance  is  known  exactly,  then  (1J) 
and  (12)  describe  rigorously  the  AJ*values  for  the 
changed  situation  as  well,  provided  we  use  the  altered 

r,*s. 

0.  The  Composite  Guide:  A  Multilayer  System 

For  a  composite  guide,  the  modes  we  have  just  dis¬ 
cussed  do  not  satisfy  the  boundary  conditions  and  can¬ 
not  propagate.  However,  they  can  be  combined  in  such 
a  way  to  produce  modes  of  the  system,  and  for  this  rea¬ 
son  wc  shall  call  them  partial  modes  to  distinguish  them 
from  the  true  system  modes.  Imagine  that  the  composite 
guide  is  formed  by  bringing  two  isolated  guides  to¬ 
gether.  I  f  these  la\  ers  are  very  fa“  apart,  then  the  modes 
of  the  combined  system  are  effectively  the  partial  inodes 
of  each  element.  However,  as  the  separation  diminishes, 
the  two  guides  interact,  and  the  modes  of  the  system 
arc  perturbed  linear  combinations  of  the  partial  modes 
belonging  to  each  layer  (cf.  Fig.  2).  Owing  to  the  sym¬ 
metry  of  the  system,  the  multilayer  modes  must  be 
either  even  or  odd  about  the  plane  of  symmetry.  For  the 
odd  modes  u(y)  must  vanish  along  the  plane  of  sym¬ 
metry,  and  for  the  even  modes  the  normal  derivative 
du/dy  must  vanish  there.  If  the  guides  are  coupled,  the 
partial  modes  in  either  layer  interact  with  one  another 
to  meet  these  conditions.  The  partial  modes  on  either 
layer  sense  the  other  and  respond  by  warping  their 
characteristics  and  their  AT- value  so  as  to  meet  the  con¬ 
dition  of  even  or  odd  symmetry.  This  can  happen  in 
two  distinct  ways,  and  since  the  reaction  would  be  dif¬ 
ferent  for  cither  parity,  the  A- value  for  e„ch  symmetry 
is  distinct.  This  discrepancy  vanishes  as  kid—*»  and 
increases  as  ktd~* 0.  This  variation  which  depends  upon 
both  »i  and  d,  is  what  we  wish  to  exploit  once  we  have 
found  its  explicit  behavior. 

The  two  layers  sense  one  another  through  the  ex¬ 
ponential  term 

exp  [±(iV»  -«,»)•'’*./]  (14) 

in  the  »t  region.  If  the  guides  are  infinitely  remote,  then 
one  can  only  use  that  exponential  for  which  ‘he  sign  is 


r,  -  y/m  -  «,*.  (13) 

For  the  TE  case,  each  I\  represents  the  normal  wall 
admittance  looking  out  of  the  guide  at  medium  »<.  For 
the  TM  polarization,  the  functions  IY*/*!*  and  T \j n»* 
are  the  corresponding  normal  wall  impedances.  Note 
th'.t  we  anticipate  the  perturbation  of  IN’*  by  affixing 
the  superscript  x.  In  (13),  both  Ti*  and  Tj.  have  the 
values  shown  because  these  are  the  required  admit¬ 
tances  needed  to  match  a  standing  wave  within  the  di¬ 
electric  slab  to  a  decaying  exponential  field  outside.  If, 

*J.  Kane  and  H.  O-terberg,  ‘Optical  characteristic*  of  planar 
guided  mode*,"  /.  Opt.  Sac.  Am.,  vof.  $4,  pp.  347-352;  March  1964. 


so  chosen  that  the  field  decays  exponentially  away  from 
the  guide.  If,  however,  the  two  guides  are  a  finite  dis¬ 
tance  apart,  then  there  will  be  a  reaction  term  and  one 
must  use  the  exponential  with  the  other  sign  to  describe 
the  coupling.  By  determining  the  field  that  must  exist 
between  the  guide  wall  and  the  plane  of  symmetry,  we 
can  find  the  eigenvalues  for  the  even  and  odd  modes  in 
the  composite  structure. 

Witn  proper  normalization,  the  expression 

exp  [+(.V*  -  y  <  0  (15) 

represents  the  unperturbed  decaying  field  from  the 
lower  wall  of  the  upper  guide  y ->  +d/2.  In  order  for  the 
total  frld  to  vanish  at  the  plane  of  symmetry  at  y*»Q, 
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we  need  to  add  the  correction  term 

-A  exp  [  —  (.V*  -  -  *)]  (16) 

to  field  (15)  whcic  the  constant  A  is 

A  -  espH-V*  “  (17) 

This  demand  of  an  odd  mode  at  the  plane  of  symmetry 
y«0  is  reflected  by  requiring  an  admittance 

1  T  ^  *_  __  .»  fit  t/ft  mill  / 1  n  \ 


r,*{rf)  -  r," 


1  -  A' 


r,«  coth  l!M(.v*  -  «i*)m]  (18) 


back  at  the  layer  boundary  y-d/2  where  IY*  i*  the 
impedance  for  *.  Likewise,  by  a  similar  calculation 
for  even  modes,  the  reflected  admittance  back  at  the 
w  alt  is  found  to  be 

IYGO  -  rr-^-4  -  r.-tanh  [}W(.V*  -  *,*)”«].  09) 

1  d-  A 

From  these  expression  we  note  the  inequalities 

r,*(<fl  >  r,-  >  IYOO  (20) 

for  any  finite  value  of  d,  and  the  relation 

iy(rf)  r,*(<o  -  (r,-)‘.  (21) 

With  a  knowledge  of  IY  and  IY,  the  eigenvalues  for  the 
multila)crcd  system  are  given  by  (11)  and  (12)  for  the 
even  and  odd  modes  provided  that  expressions  (18)  and 
(19),  respectively,  arc  substituted  for  IY*.  A  solution  of 
these  equations  then  yields  the  corresponding  iV-valucs 
for  the  even  and  odd  modes.  The  solution  of  these  im¬ 
plicit  equations  is  straightforward,  but  a  computational 
nuisance.  However,  it  is  possible  to  draw  some  conclu¬ 
sions  without  any  effort. 

An  increase  in  I\*  is  known  to  have  the  effect  of  low¬ 
ering  the  iV-value  of  the  mode.  Likewise,  a  decrease  in 
IY*  increases  the  /V- value.  From  the  inequality  (20)  we 
infer  that  the  odd  mode  has  a  higher  ,V-valuc  than  the 
corresponding  mode  of  the  isolated  guide  while  the  even 
mode  has  a  lower  one.  This  is  to  be  expected  on  physical 
grounds.  The  parameter  .Vean  be  thought  of  as  an  effec¬ 
tive  index  of  refraction  for  the  mode,  and  described  as 
some  weighted  average  of  Hi,  Hi.  and  w».  The  weight 
factors  would  be  proportional  to  the  relative  field  con¬ 
centration  in  any  one  region.  The  odd  mode  has  a  node 
in  the  «i  region,  whereas  the  even  mode  has  a  loop 
there.  As  a  result,  a  greater  portion  of  the  field  of  the 
even  mode  is  in  the  ~j  region  which  is  a  relatively  rare 
medium  and,  thus,  its  .V- value  is  lowered,  and  vice  versa 
for  the  odd  mode. 

If  changes  in  the  displacement  d  are  being  produced 
by  hydrostatic  stresses  which  produce  volume  deforma¬ 
tions.  then  the  dielectric  constants  *«<  will  be  altered. 
The  perturbed  dielectric  constants  */  can  be  found  by 
use  of  the  Lorentz-Lorenz  formula*  from  the  expression 


•  M,  Born  and  E.  Wolf.  Principles  nf  Opikt.  London:  Perga mon. 
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(**  +  2)  +  20.,’  -  \)V/V 

».  *  , _. — - - - —  (22) 

(».’  +  2)  \)V/V 

where  I'  and  V*  arc  the  original  and  jierturlxal  volumes, 
respectively.  The  effect  of  the  altered  dielectric  constant 
on  N *  and  N’  can  lie  included  by  substituting  »/  for 
Hi  in  the  characteristic  equations.  For  simplicity  we 
shall  neglect  the  perturbations  of  the  dielectric  constant, 
for  this  variation  affects  only  the  numerical  details  but 
not  the  basic  theory,  and  its  influence  can  be  included  at 
a  later  stage.  > 

111.  Thk  Hkat  Phenomenon 

A.  Theory 

In  coupled  waveguides  there  arc  paired  modes,  one 
even  and  the  other  odd,  which  are  continuous  perturba¬ 
tions  of  the  partial  modes  for  variations  in  d.  If  the 
separation  between  the  two  layers  is  very  large,  then 
the  distinction  between  the  A'-valucs  ol  the  even  and 
odd  modes  is  insignificant  owing  to  the  cx|K>ncntia! 
decay  of  A'  with  increasing  d  [cf.  (17)].  However,  as 
d  diminishes,  there  will  be  a  modest  but  significant  dis¬ 
crepancy  between  the  iV-valucs  for  the  inodes  of  cither 
symmetry.  We  shall  call  this  situation  loose  coupling  to 
distinguish  it  from  light  coupling  wherein  the  discrep¬ 
ancy  between  '.he  .V-valucs  can  be  unbounded.  For 
tight  coupling,  the  field  variation  of  the  modes  of  the 
multilayer  will  be  significantly  different  from  the  par¬ 
tial  modes  of  the  single  layers.  On  the  other  hand,  when 
the  coupling  is  loose,  the  multilayer  modes  are  essen¬ 
tially  linear  combinations  of  the  partial  modes,  and  this 
gives  rise  to  sonic  interesting  and  useful  effects. 

Let  us  denote  the  .V-valuc  of  the  even  mode  as  Ar*  and 
use  N*  for  the  odd  mode.  In  addition,  we  set 

A.V  =  .V*  -  .V.  (23) 

In  Fig.  2  we  have  sketched  the  formation  of  multilayer 
modes  as  a  combination  of  the  partial  modes.  If  we 
designate  the  y-variation  of  the  composite  mode  in. 
upper  ease  type  as  M.\’‘(y)  and  use  lower  case  letters 
m(y)  for  the  partial  modes  we  have  the  approximate 
representation 

Af.v'O')  *  ms[y  -  \{d  +  A)]  +  w.v[y  +  !(</  +  A)],  (24) 
My’iy)  "  m.vfy  -  \(d  +  /»)]  —  ro.vfy  +](</+  A)].  (25) 

These  equations  arc  valid  for  loose  coupling,  for  then 
the  partial  modes  need  not  be  significantly  perturbed 
to  form  a  true  system  mode.  Of  course,  we  can  rewrite 
these  equations  as 

*.v  [y  —  \(d  +  A)]  =  |  (3/V(y)  +  Mx*(y)  i  (26) 
M.vfy  +  Ud  +  A)]  -  \ [ifs’(y)  -  Mx‘(y)]  (27) 

which  might  seem  to  be  a  peculiar  way  of  doing  tilings, 
but  turns  out  to  very  useful.  These  last  two  equa¬ 
tions  show  that  we  can  think  of  the  field  on  the  upper 
layer  as  a  superposition  of  ^F,^,  and  My*,  and  a  field  on 
the  lower  layer  as  their  difference.  For  a  coupled  system 
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Fig.  2.  The  formation  of  the  system  modes  in  terms  of  the  partial 
modes  belonging  m  the  individual  layers. 


only  Ms*  or  ,\/,v*  are  the  possible  y- variations  of 
propagating  modes.  Any  other  field  must  be  represented 
as'  a  superposition  of  the  allowable  Ms'-'.  Supjiosc  we 
excite  the  multilayer  system  as  shown  in  Fig.  3.  That 
is,  we  illuminate  only  the  upper  layer  at  the  aperture 
plane  *<*0,  and  suppose  that  this  is  done  in  such  a  way 
that  m.\(y—d/ 2)  is  excited  in  the  vicinity  of  x  -  0.  How- 
docs  this  initial  disturbance  propagate  down  the  line 
since  mx(y—d/ 2)  is  not  the  cross  section  of  any  mode 
of  the  multilayer,  but  only  a  piece  of  one?  Equation  (26) 
furnishes  the  answer;  the  combination 
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Fig.  •>.  The  beat  phenomenon.  If  light  U  injected  into  on*  of  a  pair 
of  coupled  waveguide*,  it  will  onciitale  between  the  two  guide* 
with  a  beat  wavelength  A.  A*  the  guide*  come  closer,  the  beat 
wavelength  decrease*  and  the  output  itate  will  cycle  a*  A*.  B, 
A~,  C,  A ♦. 


the  coupling  is  loose,  if  the  beat  wavelength  is  large 
compared  to  a  wavelength.  Before  we  make  some  ex* 
plicit  calculations  describing  the  dependence  of  A  upon 
Bj  and  dt  let  us  see  how-  we  might  take  advantage  of  this 
effect. 

B.  Application 

Light  enters  the  system  as  shown  in  Fig.  3  so  that 
only  the  upper  layer  is  excited  at  x«0.  After  a  distance 
*  "  L  equal  to  an  integral  number  n  of  beat  wavelengths 

L  -  «A,  (30) 

a  photocell  samples  the  field  in  the  lower  layer.  In  the 
quiescent  state  there  is  null  output  from  the  photocell. 
Suppose  now  that  d  varies  in  such  a  manner  that  A.V 
changes  by  an  amount  sufficient  to  make  the  new  beat 
w-avciength  A'  have  the  value 


A '-(L±  *)/*. 


(3D 


»(*,  y)  -  i[3Lv'(y)  +  *<“iV'3/.v*(y)]«*»*'‘  (28) 


solves  this  initial  value  problem  for  it  is  equal  to 
m.\(y—d/2)  at  x=0  and  represents  the  transmission  of 
the  disturbance  J)(x,  y)  for  all  x,  since  3D(x,  y)  is  com¬ 
prised  of  sj-stem  modes.  That  is,  v.e  should  think  of  the 
apparatus  as  being  so  arranged  that  two  system  modes 
are  excited  at  y**0.  However,  owing  to  the  factor 
e*tiXt  ;n  (2g),  the  relative  phase  between  .1/*'  and 
Ms*  is  being  shifted  as  the  wave  travels  down  the  line. 
Indeed,  at  a  critical  distance  x~  A/2  such  that 


A/2 


to  A  A 


A 


(29) 


there  is  a  change  in  sign  in  the  coefficient  of  J/.v\  From 
(27)  we  sec  that  M.\')  represents  mx(y+d/2), 

a  field  in  but  the  lower  guide.  In  other  words,  the  en¬ 
ergy  has  been  transferred  from  on?  layer  to  the  other  at 
x-A/2.  A  repetition  of  the  argument  shows  that  when 
x  ~A,  the  y-variation  of  the  field  is  restored  to  its  initial 
state.  But  then  the  process  repeats  itself  anew,  or  the 
energy  oscillates  from  Iyer  to  layer.  A  cycle  is  com¬ 
pleted  when  x  «A,  and  for  this  reason  A  is  known  as  the 
beat  wavelength.  Since  A  depends  upon  AAr  which  in  turn 
.depends  upon  d,  measurements  of  A  can  be  used  to 
determine  changes  in  d.  Previously,  have  introduced 
the  notions  of  tight  and  loose  coupling  in  a  qualitative 
manner.  In  terms  of  the  beat  wavelength  we  say  that 


If  this  happet  3,  then  the  photocell  will  have  maximum 
reading.  Clearly,  the  photocell  output  will  be  a  sinu¬ 
soidal  response  to  fluctuations  in  d,  and  it  is  these  varia¬ 
tions  in  output  that  can  be  used  for  strain  gauge  in¬ 
strumentation.  There  is  a  continuous  exchange  of  in¬ 
tensity  between  the  two  layers  as  d  varies.  Suppose  we 
neglect  this  smooth  variation  and  classify  the  output 
as  being  in  one  of  three  states:  A,  the  intensity  equally 
divided  between  the  layers,  B,  all  the  light  in  the  top 
layer,  and  C,  all  the  light  in  the  bottom  layer.  The 
output  will  cycle  between  these  three  states  whenever 
the  relative  phase  between  M'  and  M 4  changes  by  r/4 
at  the  output  end.  Each  change  in  state  can  be  counted  I 
as  a  unit  transition.  A  count  of  the  cycling  between 
these  states  would  provide  a  digital  readout.  What  we 
need  is  some  specific  estimate  of  the  sensitivity  of  the 
proposed  instrument,  that  is  the  variation  of  AAr  as  a 
function  of  d. 

IV.  Solution  of  the  Characteristic  Equation 
A.  Procedure 

The  characteristic  equation  (II)  is  capable  of  yielding 
complete  information  describing  .V'  and  Ar*  once  the  per¬ 
turbed  values  IV(rf)  and  IY(d),  given  by  (18)  and  (19), 
are  substituted  for  Pi-.  Unfortunately,  the  resulting 
equation  is  an  implicit  transcendental  relation  for  S' 
and  A£  so  that  it  is  not  possible  to  display  any  simple 
formulas  for  their  determination.  However,  (It)  can 


IEKK  TRANSACTIONS  ON  GEOSCIENCE  Kl.l  C‘  l  KOMCS 


JCNE 


easily  be  solved  b\  numerical  techniques  and  the  re¬ 
sults  displaced  in  graphical  form. 

The  procedure  used  is  to  rewrite  (II)  in  the  form 


//  **  k,h 
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which  is  an  explicit  relation  for  the  determination  of 
the  optical  thickness  of  each  component  guide. 
No  matter  how  the  separation  D^ked  varies,  any 
thickness  AT  we  compute  must  have  this  reference  value.- 
We  can  pick  some  value  X,.  in  the  range 
and  solve  for  that  value  of  H  for  which  an  unperturbed 
layer  (£)=*«)  can  support  a  mode  which  propagates 
with  this  .V,  value.  If  D  is  finite,  this  N,  value  will  split 
into  two  values  X*  and  <V*.  each  of  which  depends  upon 
D.  We  know  that  m<X,<X,<N*<nt,  so  we  can 
guess  at  upper  and  lower  bounds  for  N'  *  given  any 
value  of  U,  it  we  substitute  these  bounas  lor  A—  into 
(32),  we  calculate  values  for  tt  which  will  be  cither  too 
large  <n  too  mull  as  cotupwcd  to  the  nlnuirt  value 
previously  computed  for  D  —  •.  The  discrepancy  be¬ 
tween  the  calculated  H  and  the  reference  indicates  how- 
poor  the  assumed  upper  and  lower  bounds  for  N'-*  arc. 
By  using  backward  interpolation,  these  bounds  can  be 
refined  to  yield  .V"  to  any  desired  degree  of  accuracy. 
In  practice,  this  procedure  it  eaSilx  programmed,  takes 
very  little  machine  time,  and  has  the  advantage  of  yield¬ 
ing  .V*  *  rigvrojtily  without  *ny  awumptiont  concerning 
strong  or  weak  coupling:  that  is,  its  validity  is  inde¬ 
pendent  of  the  distance  D. 


B.  Sample  Results 

Typical  results  are  illustrated  in  Fig.  4.  In  this 
illustration  and  »**=2.0.  For  convenience, 

the  abscissa  is  laid  off  in  units  of  T ■  1/D,  the  reciprocal 
of  D.  The  curves  show  the  separation  of  A'*  and  N‘  from 
initial  values  AT,- 1.3,  1.5,  and  1.7;  the  corresponding 
layer  thicknesses  H  are  marked  on  the  graphs  near  the 
tuning  points  ,V„  For  a  given  «V„  the  curves  N‘  and 
X*  arc  almost  indistinguishable  until  T  gets  sufficiently 
large.  In  addition,  the  separation  between  N‘  and  X* 
decreases  as  X,  increases;  this  is  to  be  expected  on 
physical  grounds  since  an  increase  in  N,  means  that  a 
greater  fraction  of  energy  propagates  in  the  Mj-dielcctric ; 
hence,  there  is  less  field  outside  which  can  sample  the 
presence  of  the  other  layer.  In  other  words,  for  a  fixed 
separation  D.  the  coupling  will  be  according  as 

X,  or  H  is 

C.  Siinificanee 

Suppose  wc  consider  the  variation  of  X*  and  A*  about 
the  value  D«X0,  as  shown  iu  Fig.  4  for  AT,- 1.3.  The 
values  of  A"*  and  .V*  arc  given  in  Tabic  I  to  five  signifi¬ 
cant  figures.  To  interpret  these  numbers,  let  us  suppose 
that  we  have  two  such  fibers  a  w  avelength  apart;  that  is, 


let  D-X»  l>c  the  quiescent  state.  If  the  system  is  one 
meter  long,  and  one  of  the  libers  is  illuminated  b\  light 
of  wavelength  691,18  mp,  then  a  lie.it  wavelength  A  is 
50  micromctcis  long  so  that  there  will  lie  2X10‘  !ieat 
wavelengths  along  the  axis  of  this  one-meter  device. 
Since  the  length  of  the  ss-ivin  is  an  integral  numlK-r  of 
beat  wavelengths,  the  exit  illumination  will  be  in  the 
same  layer  as  the  enhance  illumination.  Suppose  now 
the  distance  D  between  the  fillers  increases  from 
1.000X,  to  D  -  I.0526X*.  a  strain  of  5.26X10-’.  't  he 
fractional  change  in  A  is  8568/7236  so  that  there  are 
now  1.689X10*  beat  wavelengths  in  a  meter.  This 
represents  a  change  of  3.111X10*  beat  wavelengths; 
this  value  is  the  number  of  times  the  exit  illumination 
has  been  transferred  from  one  fiber  to  the  other  and 
back  to  the  original  layer.  Let  us  set  one  bit  to  be  equal 
to  a  quarter  change  of  state  as  the  limit  of  sensitivity. 
It  takesa  5. 26X  10~*  strain  to  produce  1.24X10*  transi¬ 
tions;  hence,  one  transition  corresponds  to  a  strain  of 
4.2 XI0-4.  In  other  words,  for  the  example  cited,  a 
system  one  meter  long  should  tic  able  to  resolve  a  strain 
of  the  order  of  10"*.  A  system  ten  meters  long  would  be 
able  to  asilw  j  sf rain  of  tht  order  I  tr*,  rtr. 

Several  factors  control  the  sensitivity:  these  include 
1)  the  optical  length  of  the  system,  2)  the  degree  of 
coupling  betw  een  the  tw  o  layers,  for  the  closer  they  are 
the  more  sensitive  is  AN  to  changes  iu  d,  and  3)  the 
value  of  X,  which  locates  the  position  on  the  dispersion 
curve  about  which  the  instrument  operates.  The  tensi- 
tivitv  increases  in  direct  proportion  to  the  length  of  the 
mstHiinf  nt.  I  lower* ,  the  variation  rA  sensitivity  with 
respect  to  2)  or  3)  is  more  complex  and  we  present  sam¬ 
ple  results  in  Table  It.  If  L  is  the  length  of  the  instru¬ 
ment  in  wavelengths,  wc  can  specify  the  smallest  strain 
3  that  it  can  resolve  as 


S 


&D/D  A 
C  AA/A  L 


(33) 


The  constant  c  equals  J  if  wc  choose  a  unit  transition  as 
being  one  quarter  of  a  complete  cycle  at  the  exit  end  as 
the  threshold  of  observation.  Table  II  illus, rates  the 
variation  of  the  sensitivity  5  as  the  coupling  and  X, 
var\.  In  Table  II  we  have  assumed,  for  example,  that 
X«  <=691.18  m p  and  that  the  length  of  the  interferometer 
was  ten  meters.  Note  that  the  strain  sensitivity  is  of  the 
order  of  10~\  if  we  decrease  D  from  Xo  to  JX0.  Of  course, 
an  additional  order  of  magnitude  or  two  of  strain  sensi¬ 
tivity  can  be  gained  at  the  exit  cud  if  the  luxury  of  digi¬ 
tal  readout  is  sacrificed  and  fractional  phase  differen¬ 
tials  between  the  two  layers  are  measured,  or  if  the 
length  l.  is  increased. 

One  disadvantage  is  apparent:  at  the  output  end 
there  is  no  simple  way  of  determining  whether  an  initial 
motion  which  produces  a  one-bit  change  co  espouds  to 
an  increase  or  decrease  iu  strain.  The  problem  arises 
from  the  fact  that  the  output  cycles  between  three 
states.  A,  B.  .and  C,  by  following  the  sequence 
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Fig.  4.  The  variation  of  A'  and  A'  with  respect  to  Tm  i/D.  In  this 
illustration  »u  — hi— 1.2  and  im  — 2.0. 


TABLE  I 


NfMMi'CAL  V,\t  rrs  MIR  A'*,  A’*.  AND  tilt'  Bui  WaVCLKNOTI!  a 
rok  tin:  Cask  ri  — 1.2,  m  —  2.0,  A«  — 1  .1,  and  //  — 0.4I8JA, 


n 

1 

AT* 

A* 

i  ** 

X\ 

. 

1 

X* 

I/A. 

Dimeri 

*ionlc« 

x# 

x# 

1.2500 

0.8000 

1.297 

i  m 

160  9 

3  5.5 

1.1765 

0.8S00 

1.296 

1.304 

177  4 

24.1 

l.llll 

0.9000 

1  295 

1.304 

10 1  3 

17.62 

1.0526 

0.9500 

1  294 

1.305 

85.68 

13.34 

! ,0000 

1.0000 

1.292 

1.306 

72.34 

10.30 

0.9524 

1.0500 

1.291 

1.307 

62  04 

8  13 

0.9091 

1.1000 

1  290 

1.308 

33.91 

6.53 

0.8696 

1.1500 

1.288 

1.509 

47.38 

5.31 

0.8333 

1.2000 

1.286 

1.310 

42.07 

TABLE  II 

The  Variation  of  the  Bfat  Wavelength  and  Strain  Sensitivity  as  a  Function  or  the  Oftical  Thickness  or  Each  Component 
Guide  and  its  Spacing.  The  Figures  Assume  a  Fiber  Length  or  Ten  Meters 


D(x.) 

AT.-I.SO.W-OWa, 

JV.-1.30, //-O  4I8A, 

1 

*- 

O 

W-0.655A, 

A  (A.) 

ixio-' 

A  (A.) 

IXIO" 

A  (A.) 

4X10" 

2.0 

1704 

458 

45400 

3360 

1450,0 

83000 

1.8 

909 

274 

7270 

1510 

50400 

8290 

1.6 

485 

164 

2930 

687 

16390 

3040 

1.4 

258 

99.7 

1190 

317 

5230 

1120 

1.2 

137 

61.4 

480 

150 

1690 

419 

1.0 

72.3 

38.4 

194 

72.4 

544 

163 

0.9 

52.3 

30.5 

123 

50.9 

309 

16J 

0.8 

37.7 

24. J 

77.8 

36.1 

175 

65.5 

0.7 

27.0 

19.5 

49.2 

25.9 

100 

42.4 

0.6 

19.2 

15.8 

30.9 

18.8 

56.4 

27.9 

0.5 

13.5 

13.2 

19.5 

13.7 

31.8 

18.7 

0.4 

9.51 

11.9 

11.9 

10.2 

17.8 

12.9 

0.3 

6.91 

8.3 

7.19 

8.07 

9.81 

9.20 

A,  B,  A,  C,  A,  B,  A,  C,  ■  •  ■  .  (34) 

Given  this  sequence,  we  have  no  way  of  knowing 
whether  an  initial  transition  front  say  A  to  B  cor¬ 
responds  to  a  positive  or  negative  strain.  Note  that  only 
the  initial  sign  is  in  doubt  if  we  include  a  record  of  the 
previous  transition  in  the  logic  of  the  output  circuitry. 
This  problem  arises  because  there  are  two  ways  by 
which  the  state  A,  equal  intensity  in  both  layers,  can 
be  reached.  There  arc  several  ways  by  which  this  degen¬ 
eracy  can  be  eliminated.  A  calibration  strain  artificially 
produced  can  eliiniiialc  subsequent  directional  ambigu¬ 
ity.  Another  niclhod  is  to  simultaneously  illuminate  the 
system  by  two  distinct  monochromatic  sources  and 
compare  th<*  system  response  to  either  color.  This  also 
has  the  potential  advantage  of  increasing  sensitivity. 


Instead  of  using  two  different  colors  an  equivalent  and 
simpler  procedure  would  be  to  excite  both  TE  and  TM 
modes  and  analyze  the  output  through  polarizers.  This 
would  also  provide  additional  resolution,  since  the  TE 
and  TM  modes  are  uncoupled  and  propagate  with  dif¬ 
ferent  phase  velocities.  Assume,  for  example,  that  these 
two  velocities  didet  by  ten  percent.  This  difference  can 
be  used  at  the  output  end  to  provide  a  vernier  measure¬ 
ment  and  increase  the  sensitivity  by  an  order  of  magni¬ 
tude.  A  third  alternative  is  to  use  a  thrcc-layci  system 
(sec  Appendix)  which  provides  an  additional  output 
that  can  be  used  for  directional  calibration. 

Perhaps  the  simplest  method  is  to  recognize  that  a 
distinction  can  be  made  in  state  A,  for  when  the  inten¬ 
sity  is  equally  divided  between  the  two  layers  there  is  a 
ir/4  phase  difference  between  the  fields  in  the  top  and 
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bottom  l.i \ cm.  As  i  refill*.  A  t’.ni  Ik*  divided  into  two 
substates:  A\  wliit  h  nuurs  when  ti  ■  phase  in  the  top 
layer  leads  the  ph.ow  in  the  Ixrttotn  layer  l>\  tr/4.  and 
A  ,  w  hicli  irises  w  lieu  the  pli.isc  in  the  top  Inver  lags  t  he 
phase  n  the  lower  liver  by  »/4.  llencc,  if  phase  dis¬ 
crimination  is  included  in  the  scnsni*,  and  this  can  lie 
done  easily,  then  the  response  would  cycle  as 

A\  B,  A  ,  C,  A\  B,  A~,  C,  A*,  •  •  •  .  (35) 

\\  i(h  ain  of  these  refnieinents,  a  transition  between  any 
two  adjareitt  states  would  leave  no  doubt  as  to  the  sign 
of  the  strain. 

V  rV.sir.N  Notks 

.4.  Fabr nation 

Geologists  arc  a»  icd  to  thinking  of  layered 

structures  .is  being  r.iti.ci  .urge.  For  simplicity,  wc  have 
assumed  that  the  r-dimension  of  our  waveguides  is  in¬ 
finite.  In  practice,  it  is  only  necessary  that  this  dimen¬ 
sion  be  large  compared  .o  a  wavelength.  At  optical  fre¬ 
quencies  this  criterion  can  be  achieved  w  ith  a  depth  of 
but  a  few  millimeters. 

On  the  other  hand,  the  problem  of  achieving  dimen¬ 
sional  tolerances  may  appear  to  be  formidable.  This 
would,  of  course,  be  true  if  the  components  of  the 
dielectric  guides  were  manufactured  separately  and 
then  assei  iblcd  in  the  microscopic  state.  However, 
fiber  optic  devices  arc  not  made  this  way;  they  are 
fashioned  from  rather  large  glass  samples.  All  of  the 
components  are  assembled  on  a  macroscopic  scale  in  a 
glass  matrix  with  an  extreme  magnification  of  scale  in 
the  cross  sectional  dimensions.  Following  this,  the  com¬ 
posite  material  is  heated  and  drawn  into  a  fine  thread 
which  shrinks  the  cross  section  by  several  orders  of 
magnitude*.  F.xccpt  for  a  scale  change,  the  cross-scc- 
tional  integrity  is  maintained  as  the  fiber  is  drawn.  In 
this  manner,  current  technology  has  no  difficulty  in 
producing  composite  glass  fibers  with  square  cores. 
These  have  been  used  in  fiber  optic  lasers  to  increase 
their  packing  density  and  hence  efficiency.  Conse¬ 
quent!).  there  should  be  no  problem  in  fashioning  not 
only  a  pair  of  waveguides,  but  a  w  hole  array  of  them  at 
once  by  drawing  a  multicomponent  glass  matrix.  For 
example,  this  allows  the  possibility  of  drawing  many 
waveguides  at  once.  In  the  configuration  of  Fig.  5,  pairs 
of  coupled  waveguides  are  at  right  angles.  By  mixing 
the  output  in  appropriate  linear  combinations  it  is  possi¬ 
ble  to  devise  a  strain  gauge  w  hich  can  discriminate  be¬ 
tween  transverse  shear  and  volume  deformations. 

A  review  of  the  analysis  will  indicate  that  tolerances 
in  dimensions  are  not  critical,  for  w  hat  is  being  measured 
is  the  average  strain  in  a  direction  normal  to  the  axis  of 
the  dielectric  waveguides.  In  fact,  this  averaging  of 
strain  information  should  also  have  the  beneficial  effect 
of  reducing  noise  level. 


Fig  5,  .*  diagram  rtf  a  gla*s  matrix  eonlainiiiK  an  arr.i>  of  wave- 

guidex  ^rior  to  drawing  I  ol lowing  the  drawing  proccvt  the  pair, 
of  fibers  can  )  ield  information  on  twro  components  of  transverse 
strain. 


H.  Directional  Response 

The  proposed  device  has  a  directional  response  for  it 
responds  essentially  to  strains  in  the  y-dimensiott  only. 
It  is  not  hard  to  see  that  its  polar  pattern  is  basically  a 
cos$  pattern  in  the  x,  y  plane  w  ith  the  major  lobe  in  the 
direction  of  the  y-axis.  The  strain  rejection  in  other 
dimensions  is  due  principally  to  the  fact  that  the  only 
response  that  can  even  be  elicited  occurs  when  a  dis¬ 
placement  produces  a  differential  between  N"  and  N‘. 
That  is,  even  though  a  strain  may  change  both  ,Vr  and 
N’,  the  output  will  respond  only  to  strains  which  alter 
iV*--.V*  and  not  to  motions  which  change  iV'  +  JV*. 

C.  Other  Configurations 

Although  the  device  wc  have  described  is  a  self-con¬ 
tained  transverse  strain  gauge,  the  ideas  inherent  in  its 
operation  can  be  used  in  other  useful  ways.  Laser  inter¬ 
ferometry  has  had  high  promise  in  measuring  small  geo¬ 
logical  strains,  but  this  anticipation  has  been  largely 
aborted  owing  to  the  unstable  character  of  atmospheric 
propagation  paths  w  hen  the  optical  paths  arc  large.  The 
need  here  is  ft:  interferometry  along  a  sheltered  optical 
circuit.  Dielectric  lines  offer  one  solution  to  this  problem 
by  providing  a  stable  propagation  path.  Figure  6  illuc- 
t rates  a  Michelson  interferometer  devised  on  the  prin¬ 
ciples  just  described  to  fashion  a  longitudinal  strain 
gauge.  Traditionally,  the  Michelson  interferometer  re¬ 
quires  a  half-silvered  mirror  and  two  mirrors  to  reflect 
the  light  at  the  end  of  two  perpendicular  arms.  Wc  can 
achieve  these  requirements  in  the  following  fashion.  A 
pair  of  coupled  fibers  at  A,  w  hich  is  («  +  «)  beat  wave¬ 
lengths  long,  functions  as  a  directional  coupler  and 
divides  a  coherent  input  into  two  branches  after  the 
fashion  of  a  haif-silvered  mirror.  The  fibers  then 
separate  into  two  long  paths  of  lengths  /t  and  It  at  right 
angles.  Instead  of  a  reflecting  surface,  the  dielectric 
lines  can  be  bent  at  B  by  180°  provided  the  radius  of 
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Fk-  ••  A  Michcl'nii  inlcrfi-rniuclcr  (used  on  fiber  optic  waveguide 
’  riuciplev  The  liber  waveguides  are  coupled  only  in  the  ertws- 
alched  regions  A  and  C.  A  beam  of  lighl  in  one  fiber  will  be 
divided  between  I  wo  waveguides  if  the  cotipler  A  is  so  adjusted 
that  il  is  (n  +  JI  beat  wavelengths  long.  After  traversing  two  arms 
of  ihe  interferometer,  the  light  is  bent  at  B  and  broughl  back 
to  a  long  coupler  C.  The  phase  difference  between  the  two  beams 
will  be  mulliplied  by  a  factor  depending  upon  the  length  of  L 
before  being  analyzed. 

curvature  r  is  greater’  than  r«,i„ 

(35) 

The  two  fibers  arc  then  coupled  together  at  C  to  provide 
a  comparison  of  the  path  lengths  h  and  /j.  If  the  length 
L  of  the  coupled  fibers  at  C  is  long,  then  amplification  of 
the  phase  difference  produced  by  variations  in  l\  —  lt 
takes  place  and,  as  described  previously,  digital  readout 
can  be  achieved  at  D. 

D.  Comparison  with  Existing  Strain  Gauges 

'rhe  proposed  instruments  have  several  unique  fea¬ 
tures  which  prevent  a  meaningful  comparison  with 
available  strain  gauges.  Existing  devices  as  a  rule  mea¬ 
sure  the  strain  between  a  pair  of  points.  For  example, 
the  Benioff  strain  seismometer  measures  the  relative 
motion  of  two  piers,  and  a  crystalline  transducer  mea¬ 
sures  the  sttain  of  two  of  its  facets.  The  transverse  strain 
gauge  we  have  described  measures  the  parallel  displace¬ 
ment  of  two  lines  rather  than  the  distance  between 
points.  For  this  reason,  the  two-fiber  system  is  more 
comparable  to  an  array  of  strain  gauges.  Therefore, 
while  the  theoretical  strain  sensitivity  we  have  calcu¬ 
lated  is  not  competitive  with  the  sensitivity  of  a  first- 
class  Benioff  installation,  its  readings  should  have  a 
lower  noise  level,  inasmuch  as  we  are  reading  parallel 
displacement  axeraged  over  a  finite  distance.  To  appre¬ 
ciate  this  point,  consider  the  parallel  displacement  re¬ 
quired  to  yield  the  figures  in  Tables  1  and  1 1  for  a  pair 
of  fibers  ten  meters  long  and  half  a  wavelength  apart 
where  X»-  100  nip.  If  we  assume  that  the  system  param¬ 
eters  arc  such  that  the  second  column  in  Table  II  is 

*  S.  E.  Miller,  "Directional  control  in  light-wave  guidance,'  Bril 
Sys.  Trek  J ..  vol.  45,  pp.  1727-1739,  july  1964. 
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applicable,  then  we  have  a  strain  sensitivity  of  14.7 
XIO  ’.  In  trims  of  displacement,  this  means  that  we 
should  lie  able  to  resolve  a  motion  of  alxjut  0  003 
angstrom,  a  distance  in  itch  smaller  than  the  diameter 
of  .lie  hydrogen  atom.  Clearly  any  point-to-point  dis¬ 
placement  of  this  magnitude  would  be  impossible  to 
measure  owing  to  quantum  and  thermal  effects.  How¬ 
ever.  there  is  no  theoretical  reason  that  limits  the  resolu¬ 
tion  of  the  average  displacement  o!  two  lines.  Indeed, 
the  transverse  displacement  sensitivity  of  the  filler  pair 
implies  that  the  device  might  have  other  geological 
uses  than  a  strain  gauge.  In  A  conversation,  Dr.  F.  I’rcss 
has  suggested  that  its  most  immediate  application 
might  lie  as  a  digital  displacement  transducer  to  mon¬ 
itor  the  small  air  gap  in  a  Benioff  seismometer. 

VI.  CONTUSIONS 

In  this  study  we  have  introduced  several  novel  appli¬ 
cations  of  filter  optics  to  the  problem  of  optical  seis- 
nioinctrv.  Although  the  particular  configurations  de¬ 
scribed  may  not  prove  to  be  competitive  with  existing 
instruments,  the  ideas  incorporated  in  their  operation 
may  prove  to  have  unexpected  geophysical  applications. 
Several  features  stand  out  and  arc  .vorthy  of  emphasis. 
First  is  the  method  of  achieving  digital  readout  by 
using  a  beat  phenomenon.  This  also  has  the  advantage 
of  operating  the  device  in  a  balanced  state,  so  that  the 
effect  of  local  perturbations  are  averaged  out  in  the 
output  stage.  As  a  result,  the  noise  level  of  the  fiber 
optic  strain  gauge  should  be  lower  than  current  devices. 
The  basic  simplicity  is  another  feature,  and  this  would 
be  reflected  in  a  relatively  modest  instrument  cost.  In 
addition,  the  design  is  such  that  the  device  could  oper¬ 
ate  in  hostile  environments  owing  to  the  inherent 
strength  and  corrosion  resistence  of  glass  fibers. 

Appendix 

The  Three-Layer  System 

A.  Analysis 

The  addition  of  a  third  layer  as  indicated  in  Fig.  7 
can  provide  an  unambiguous  three-stage  output  once 
we  understand  the  characteristics  of  mode  propagation 
in  this  system.  The  simplest  way  to  describe  the  modes 
of  the  three-layer  system  is  to  consider  them  as  linear 
combinations  of  the  modes  .1/,*,  3 f.v*  of  the  two-layer 
system  which  are  partial  modes  of  the  three-layer  con¬ 
figuration.  Let  us  designate  the  y-variation  of  the  sys¬ 
tem  modes  as  The  possible  system  modes 

9U.v(y)  can  be  thought  of  as  perturbed  linear  combina¬ 
tions  of  the  two-layer  modes  .If.,-*  and  3f.v*.  The  equi¬ 
distant  configuration  sketched  in  Fig.  7  possesses  a 
plane  of  symmetry  about  the  cenfial  layer,  hence,  the 
Sffl-modes  wil!  be  even  and  odd  combinations  of  its* 
and  The  partial  mode  3 f.v'  can  be  combined  in  two 
ways  to  form  an  311-mode,  one  even  and  the  ether  odd. 
vi*. 
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w.v"(y)  -  1  |3M>- W  +  *)l 
+  A/.v‘[y+i(rf  +  *)l| 
sn,v**o>)  «  J|A/,v'[y  -  \(d  +  *)] 

-  MSb+Ud  +  k)  jj 

and  likewise,  .Wv*  can  be  arrange  J  in  two  ways, 

v  -  Wifs’b-Kd  +  k)] 

+  Sh*b-k(d  +  h)\ 

SU.v-  «=  J  { M.v*[y  -  l(rf  +  A)]  ‘ 

-KAy-\(d  +  h)  ] 


We  can  rewrite  these  expressions  in  terms  of  the  y-varia- 
tion  of  the  primary  modes  m(y), 


sm,v"(y)  «  i{m.v[y  —  !(d  +  A)]  +  2ms(y) 
+  m.v[y  +  W  +  A)JJ 

9K.v**(y)  -  3TC.v'*(y)  -  l{«.v[y  --  W  +  A)] 

-«.v[y+K<*  +  *)]} 

3u.v"Cy)  ■■  l{«.v[y  -  W  +  *)]  ~ 2 m"(y) 

+  «.v[y+ 1(<M- *)]} 


.  (38) 


Note  that  3H**  and  M'*  are  degenerate,  that  is  311**  —  UTC**, 
cr  they  represent  the  same  state. 


T.f/O— ► 


Fig.  7.  The  geometry  of  a  three- layer  device  and  value*  of  N". 
A'*,  and  .V*  a*  a  function  of  T-l/D  for  the  case  1.2, 

*•-2.0. 


B.  The  Characteristic  Numbers 

The  eigenvalue  N  will  be  different  for  each  of  the 
modes  Olt.v",  311  .v",  3%"*,  and  we  can  designate  these 
three  values  as  .V",  A'*",  and  Vw,  rcsjiectively.  Appro  ti¬ 
nt  a  ic  characteristic  equations  describing  these  vain  :s 
arc  readily  found  by  concentrating  attention  on  the 
central  layer.  To  the  first-order,  the  three  numbers  A’**, 
A**,  and  Nm  represent  the  three  ways  by  which  the  wall 
admittances  I\  of  this  layer  can  be  perturbed.  We  can 
combine  the  possible  perturbations  to  find  the  charac¬ 
teristic  equations  v.  hich  describe  the  three-layer  system 

tan  {[  ,*  -  (.Vi0,],/’AoAj 

„  K  -  Qv^MivM  +  lYffll  TE 

»*’-  r,‘(d)-rt'(d)  ’  {  } 

in  this  expression,  the  superscripts  i,j  index  the  possible 
states  t  =  e  o  and  j=e,  o.  It  might  seem  invalid  to  calcu¬ 
late  A"  from  an  argument  involving  the  field  in  the 
central  layer  w  hich  presumably  is  null,  i  low  ever,  wc 
remind  the  reader  that  the  field  there  will  be  relatively 
small  and  not  null  since  311"  is  comprised  not  of  m.v(y) 
modes,  but  rather  of  somewhat  perturbed  versions  so 
that  complete  cancellation  of  field  never  takes  place  in 
the  central  layer.  Even  though  the  resulting  field  in  the 
layer  is  very  small,  (39)  can  be  used  to  find  A"  since 
(39)  depends  only  upon  the  logarithmic  derivative  of  the 
field  which  is  independent  of  amplitude  considerations. 

A  feature  of  the  point  of  view  just  introduced  is  that 
it  allows  us  to  readily  estimate  the  values  A"(</),  N‘°(d), 
and  N°*(d)  in  terms  of  the  preceding  analysis.  The  num¬ 
bers  Nu(d),  N°*(d)  coi  respond  to  a  tw  ofold  per¬ 

turbation  of  the  wall  admittances  of  the  central  layer. 
To  first-order,  the  two  perturbations  are  independent  so 
that  we  can  combine  the  A*(d)  and  N~(d)  values  of  the 
two-layer  system  to  obtain 

A«(<0  -  AT,  +  2[.V.  -  .¥•(</)] 
iV**(<f)  -  Ar,  +  Ar,(<f)  -  A T*(rf) 

AT"(<0  -  -V.  +  2[.Y’(rf)  -  AT.],  ■  (40) 

which  are  estimates  for  AT”,  Ar**,  and  Ar**.  In  terms  of 
this  simple  argument,  wc  can  plot  the  typical  behavior 
of  AT",  A",  and  A**  (cf.  Fig.  7)  and  note  that 

A '"(4  <  A”(<0  <  N"(d).  (41) 

C.  The  Beat  Phenomena 

If  we  invert  the  set  of  (38),  we  find 

*.v[y  -  J(rf  +  A)]  -  3Jl.v"(y)  +  23 TC.v**(y)  +  tra.v"(y) 
«,v(y)  *  3TC.v"(y)  -  3 tt.v“(y) 

w.v[y  +  \{i  +  A)]  -  DU.v**(y)  -  231t.v"(y)  +  3U,V~(>).  («) 
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The  sini.ifK-u.ee  nf  u  -e  equations  is  th.it  they  allow 
«s  to  exp  ess  illtimiuati  >n  in  any  one  layer  as  a  linear 
eombiiiatioi:  o(  three -system  moden.  Since  the  Ar  values 
o(  tliesf*  modes  are  distinct,  the  relative  phases  of  the 
imxlc  that  are  excited  vary  as  a  function  of  distance 
along  the  three-layer  combination.  Tints,  if  we  illu¬ 
minate  only  the  top  layer  at  the  entrance  the  light  will  be 
transferred  In-tween  the  layers  as  the  disturbance  propa¬ 
gates  down  the  line. 

As  an  example,  Jet  us  choose  fti-1.2,  »j“2.0,  H 
*=Q.41SX(,,  and  D~\ 9.  Fron  Table  I  and  (40),  we  can 
use  the  values 

X”  *=  1.286,  A'**  «=  1.298,  N"  -  1.312  (43) 

for  the  t’lrec-laycr  system.  If  we  inject  illumination  only 
at  the  top  layer  at  .v=0,  the  expression  describing  the 
propagation  of  the  disturbance  »(.r,  y )  is 

»(*,  y)  =  +  23TlJv*‘(j)e,  ***i*»* 

4-  3lt.v"(y)e"l,,*M  (44) 

or 

45>(*t  y)  -  (e*  «**•*  +  2e*  «•«** 

+  el",,‘M]w.v(y  —  d) 

||i.itiaii _ 2r‘  nta%“ 

+  +  d).  (45) 


The  beat  phenomenon  that  now  takes  place  is  a  good 
deal  more  complex.  In  Fig.  8  we  plot  the  relative  in¬ 
tensity  o?  illumination  in  each  layer  as  a  function  of  dis¬ 
tance  kiyX  from  ?;cro  to  lOOOr  or  !00  wavelengths.  This 
interval  returns  the  system  to  its  original  state  apart 
from  a  sign  change,  and  the  cycling  begins  anew.  Note 
the  basic  difficulty  of  the  three-layer  system  which  is 
the  apparent  apcriodicity.  This  results  from  A7**,  N**, 
and  A'**  having  a  rather  large  least-common  divisor.  For 
a  useful  instrument,  it  would  be  necessary  to  devise  a 
three-layer  configuration  such  that  N",  N“,  and  N~ 
follow  some  elementary  ratio  like  4:5:6. 


Length  {X0)-* 


The  beat  phenomenon  in  a  system  of  three  waveguides; /■, 
ft,  and  ft  are  proportional  to  the  itUtntUy  of  light  if  the  excitation 
is  in  the  top  guide  at  x-0. 
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Vector  Wave  Diffraction  at  Crustal  Discontinuities 
L  Basic  Theory,  Rayleigh  Waves  at  the  Continental  Margin 

Julius  Kane* 


Summary 

....  Jo  this  report,  the  first  of  a  series  on  the  analytic  continuation  of  wave 
functions  past  discontinuities,  we  introduce  an  elementary  procedure 
for  the  solution  of  problems  involving  the  diffraction  of  vector  fields. 

In  particular,  the  report  discusses  the  propagation  of  Rayleigh  waves 
incident  obliquely  upon  the  continental  margin.  The  crustal  layering 
on  either  side  of  the  coastline  is  modelled  mathematically  as  a  two-part 
boundary  layer  in  such  a  fashion  that  the  relevant  reflection  and  trans¬ 
mission  coefficients  emerge  as  elementary  algebraic  expressions.  The  pro¬ 
cedure  permitting  such  a  solution  to  be  found  is  to  introduce  a  diffraction 
analogue  of  the  well-known  procedure  in  electrical  engineering:  replace 
a  transmission  line  network  by  a  lumped  parameter  equivalent  in  a 
specified  frequency  or  wave  number  domain. 

I.  Introduction 

From  the  point  of  view  of  a  geophysicist,  mathematical  analysis  cannot  begin 
to  acquire  practical  importance  unless  the  solution  is  in  such  a  form  that  the 
answers  to  his  questions  can  be  promptly  given.  Unfortunately,  the  class  of  elastic 
wave  diffraction  problems  which  have  tractable  solutions  usually  describes  over¬ 
simplified  geophysical  structures.  On  the  other  hand,  when  the  theoretician  studies 
complex  problems,  too  oil  n  he  concentrates  on  the  mathematical  details  or  the 
solution  so  that  his  calculations  are  often  more  impressive  than  useful. 

One  powerful  mathematical  technique,  the  Wicncr-Hopr  method,  has  not  had 
the  application  to  geophysics  that  it  should  have.  The  difficulty  in  the  application 
of  this  method  is  that  it  usually  requires  an  involved  transcendental  kernel  to  be 
•split*  into  suitable  analytic  factors.  In  my  experience  with  problems  A  this  type, 
it  has  been  my  observation  that  very  often  this  prolix  calculation  is  of  more  mathe¬ 
matical  than  physical  interest  (1,  2),  and  that  usually  the  interesting  parts  or  the 
solution  depend  in  a  very  minor  fashion  with  respect  to  this  ’split*  function.  Indeed, 
the  necessity  for  this  function-theoretic  decomposition  can  often  be  dropped  as  a 
requirement  from  the  analysis,  and  used  as  the  basis  of  an  approximation  procedure 
to  solve  problems  beyond  the  scope  of  the  Wicner-Hopf  theory  (3).  Other  investi¬ 
gators  such  as  Koiter  (4)  and  Carrier  (5)  have  been  aware  of  the  relative  unimportance 
of  this  complex  factorization  and  have  suggested  the  use  of  simpler  substitute 

•Present  address:  Space  Science  Center,  University  of  California,  Los  Angeles.  California  U.S.A. 
(on  academic  leave  from  the  University  of  Rhode  Island,  Kingston,  Rhode  Island,  U.S.A.). 
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kernels  which  can  be  fa-.  lorizxd  by  inspection.  While  there  is  much  merit  in  this 
approach,  I  feel  that  it  is  a  poor  philosophy  to  take  the  solution  to  an  idealized 
problem  and  then  approximate  the  answer.  A  preferable  procedure  is  to  formulate 
the  problem  in  such  a  way  that  the  mathematical  difficulties  are  anticipated  and 
avoided  from  the  start.  That  is,  a  certain  amount  of  intuition  and  insight  into  the 
nature  of  the  problem  can  be  used  to  guarantee  that  the  subsequent  analysis  will 
lead  to  simple  representations. 

In  the  sequel,  I  should  like  to  give  an  example  of  this  philosophy  by  considering 
the  propagation  of  oblique  Rayleigh  waves  past  a  crustal  discontinuity  such  as  the 
‘  .nd-sca  boundary  (Fig.  1).  As  a  side  condition,  an  elementary  character  for  the 
i  ution  will  be  required.  Of  ccursc,  thete  will  be  a  need  to  make  certain  assump- 
uo.it  concerning  the  fields  to  b  calculated.  However,  we  shall  be  explicit  about 
these  modifications  and  indicate  how  arbitrary  improvements  can  be  made  in  the 
analysis  if  necessary. 


Fig.  t.  A  Rayleigh  wave  is  incident  obliquely  upon  the  land-sea  interface  at 
an  angle  The  transmitted  wave  is  refracted  by  an  angle  r.  For  purposes  of 
illustration,  the  crustal  thickness  is  greatly  exaggerated  and  is  considered  to 
be  a  two-part  boundary  layer  in  the  analysis. 


A  crucial  phase  in  this  study  will  be  the  need  to  formulate  appropriate  boundary 
conditions  to  characterize  the  crustal  layering  of  the  mantle.  In  the  simpler  problem 
of  electromagnetic  wave  propagation,  Leontovich  (6)  has  introduced  the  notion 
of  an  'impedance  boundary  condition’  to  model  the  physics  at  an  interface.  How¬ 
ever,  this  idea  while  very  useful  does  have  its  limitations.  For  one  the  solution  to 
the  associated  boundary  value  problem  is  still  rather  complex,  second  the  idea  is 
valid  only  if  the  surface  impedance  is  reasonably  independent  J  the' nature  of  the 
excitation.  Although  the  impedance  concept  can  be  generalized  to  describe  inter¬ 
laces  with  more  general  properties  (7-9),  the  intractable  nature  of  the  solution 
remains.  A  feature  of  the  subsequent  analysis  will  be  a  derivation  of  an  approximate 
boundary  condition  for  elastic  wave  diffraction  whose  merit  is  that  its  accuracy 
can  be  arbitrarily  great  without  sacrificing  simplicity  in  the  form  of  the  solution. 
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2.  Basic  equations  .  ... 

In  the  absence  of  body  forces,  small  displacements  s (».  v.  w)  of  an  elastic  solid 

characterized  by  the  Lame  parameters  A,  P,  and  density  p  can  be  derived  from  a 

scalar  potential  <b  and  a  vector  potential  Y(Y„0,H',) 

s  -  grad «b+ curin'^,, 0,H»j)  (0 

(Three  components  in  the  vector  potential  arc  redundant  for  we  can  express  say 
u#  «=  rtY,  Y\)  so  that  without  any  loss  of  generality  we  O'.!  set  Y*  »  0.)  ror 
monochromatic  vibrations,  we  can  oppress  a  time  factor  -  and  it  can  be  shown 
that  4>  and  the  components  Y,  satisfy  the  reduced  wave  equations 

(Vi+kf)G>  *=  0,  k*  =  <oxpl(X+2p) 

(V2  +  k,2)'Fi  =  0,  k,2  =  «V/P- 

Simultaneously  we  should  like  a  solution  for  the  Rayleigh  wave  incident  from 
either  side  of  the  crustal  discontinuity.  For  this  reason  let  J*  ,hc  Pja*« 

velocity  of  the  Rayleigh  wave  on  the  left  and  right  sides  or  the  boundary,  and  let 
kL  m  alp.,  and  *,  -  afvn  designate  the  corresponding  wave  numbers.  Also, £ 
and  r,  represent  the  required  shear/compressional  ratios  for  either  side.  In  the 
notation  Fig.  1,  the  variation  of  the  incident  Rayleigh  wave  in  the  mantle  will 

have  the  form  (10) 


where 


l  Yu..  “  (««  e-  0.  cos  ©)tfy 


{  s  xcos©+zsin  0, 


a  =  (kl-k?)*,  b  =  (kl-W,  (5) 

and  for  the  special  case  of  a  free  elastic  half  space  the  coefficient  T,  would  be  given 

by 

[l-(unft,)*l*  (6) 

The  z-variation  of  all  fields  will  be  of  the  form  i.e. 

(7) 

W*'  y>  *)  -  w*.  y*'  1  *  3 

so  that  «*x,y)  and  the  satisfy  the  reduced  two-dimensional  wave  equations 

p7  =  ti-kl  sin2©, 

dx*  dr  .  ’  (8) 

0.  q1  *  fcj— fcjsin2  ©. 

Note  that  for  angles  of  incidence  such  that  sinO  >  vjv,  both  p2  and  q1  will  be 
negative. 
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A  completely  general  solution  of  equations  (8)  can  be  written  down  at  once 
by  separation  of  variables  as 


y)  =  —■  |/l(at)  :xp[»auc  -  i(p2  -  aJ)»y)da 


y)m2ni  J  S‘(3>  cxpliajc “ <9* ~  **>**)<*■>  i  »  1, 3 


(9) 


where  /1(a)  and  B&a)  are  kernels  to  be  found,  and  the  contour  V  follows  the  real 
axis  except  for  the  standard  pole  and  branch  point  deformation  (Fig.  2). 


Fig.  2.  The  complex  a-plnnc.  The  contour  follows  the  real  axis  except  for 
the  indicated  pole  and  branch  point  deformations.  The  pole  at  I  does  not 
play  any  role  in  the  analysis  and  is  not  shown. 

3.  The  boundary  conditions 

Either  the  oceanic  or  the  continental  crustal  structure  can  be  considered  to  be 
a  shallow  transition  layer  that  continues  the  field  within  the  mantle  to  the  free 
boundary  at  the  Earth’s  surface.  The  details  of  the  transition  become  more  important 
for  the  shorter  period  Rayleigh  waves.  However,  the  group  velocity  dispersion 
curves  for  Rayleigh  waves  travelling  oceanic  or  continental  paths  intersect  at 
about  17  s.  This  implies  that  for  waves  of  period  greater  than  10  s  we  can  treat 
such  a  crustal  structure  as  a  two-part  boundary  layer  whose  coefficients  have 
different  values  on  the  land  and  sea  sides  of  the  coastline.  For  this  purpose,  we 
need  some  way  of  relating  the  free  surface  conditions  to  the  crust-mantle  interface. 

A.  Exact  fomUiC*:on.  At  a  free  surface,  the  normal  and  tangential  stresses  must 
vanish-or 

(diu  <)*w\  to 


du  dv  . 

‘ 

dw  dv  „ 

dy  di 


(10) 
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which  in  terms  of  <p  and  tfq  become 


d*(p 


2-f+i*<tSin0,fi  + 


,  r’Vj  t’*Vj 


"  dxdy  *  dx  ,V  dxJ  ~ 


(ID 


By  means  of  the  Haskell-Thomson  (11,12)  malrix  method,  these  equations  at 
the  surface  can  be  transformed  to  an  equivalent  set  which  refer  the  boundary 
conditions  a'  the  free  surface  to  the  crust-mantle  interface.  In  this  fashion  the 
details  of  the  crustal  layering  on  either  side  of  the  boundary  will  be  introduced  into 
the  analysis.  However,  the  boundary  conations  so  obtained  will  be  no  simpler 
than  the  set  (1 IX  and  in  fact  a  good  deal  n;w«c  complex  Although  the  programme 
we  have  just  described  would  lead  to  a  soluble  Wicncr-Hopf  problem,  the  com¬ 
plexities  of  the  analysis  would  render  this  procedure  rather  prolix  for  the  information 
desired. 


B.  The  simplified  formulation.  The  important  features  of  the  exact  foimulath"- 
described  in  part  A  can  be  inr  orated  in  a  muc  »  i.nplcr  r  ,  >ach.  The  condi¬ 
tions  we  need  arc  those  that  characterize  the  boundary  1  ycr  h  ,uch  a  fashion  that 
the  geometric  acoustics  poles  mirror  the  crustal  layering  correctly.  We  can  assume 
that  we  know  the  dispersion  relations  for  cither  side  or  the  boundary,  this  i' 
equivalent  to  making  use  of  the  Haskell-Thomson  matrix  analysis  in  an  indirect 
fashion. 

The  dispersion  relations  define  the  vK  L  and  Vf  L  of  a  Rayleigh  wave  un."  Ay. 
If  we  refer  to  the  explicit  representation  (3)  we  see  that  at  any  characteristic  depth 
say  y  «  0  the  Rayleigh  wave  will  satisfy  the  set  of  equations 


.  ^  A  d'l\  <?4'3 

km  I  r a  lO  —  - _ — 

*.L  K.L 


dZ  Ijj. 


g+^iVj=.o. 


(12) 


These  equations  arc  analogues  of  the  telegrapher’s  equations  in  electromagnetic 
theory,  and  arc  used  to  describe  the  propagation  characteristics  of  Rayleigh  waves 
in  a  layered  half  space.  They  describe  the  surface  loading  of  the  mantle  in  an  exact 
fashion  for  Rayleigh  wave  propagation.  That  is,  they  define  the  correct  velocity 
pHit,  and  the  shear/comprcssional  ratio  ±  T F,L  together  with  its  proper  change  of 
sign  when  the  direction  of  propagation  is  reversed.  The  surface  parameter  r,,it 
is  depth  dependent  and  transforms  according  to  the  taw 
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(13) 
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in  velocity  (Fig.  3). 


M  ...  r>  jJ _ _ D„v — . 

SjSS  ,  ;N  ;R«f«r«w*  pUn# 


tr,„  •>  An  eusljc  soliJ  wilh  a  small  step  discontinuity  can  be  thought  of  as 
SL,  v.h*  of  F  or  ..hr  ,*  of  .h  «» 

analysis  is  referred  to  the  reference  plane  indicated. 

or  course  there  arc  any  number  of  other  simplified  boundary  conditions  equa¬ 
tion that  the  Rayleigh  wave  will  satisfy,  but  the  equations  m  the  set  (12,  have 

several  important  advantages.  ( 

fa,  They  ate  simple  and  uncomplicated. 

(b)  They  are  isotropic  in  that  they  do  not  depend  upo.  tc  y 

direction  of  propagation  in  the  xr-planc.  havc  ^ 

/ct  Thcv  avoid  the  use  of  the  transverse  operator  c/cy  ....  -j 

effect  of  introducing  radicals  in  the  transform  plane.  This  j 

ance  of  the  irrational  factors  will  enormously  simplify  the  subseque 

function-theoretical  analysis.  ..  lhev  ;ntroduce 

«  t,  ray, 
imposed  will  have  the  same  energy  as  the  incident  plane  wa  . 
c.  D„To,  For  .hose  Wf,  »"t‘SuX  CST-SS 

ESd"iT.rix“  mrthX  a  vector  aLloE«e  ot  .ransforming  an  impedance  throng 

line  representation  or  the  mantle.  As  is  welt  known,  the  impedance  characteris 
nf  such  networks  are  involved  transcendental  functions  of  irequcncy.  Within  ^y 

however,  a  transmission  line  network  c-tejjjjd 

by  an  equivalent  network  consisting  of  lumped  parameter  elemen  h 

2‘ much  simpler  impedance  functions  which  are  rational  funct.cn*  of  he 

frequency.  In  diffraction  problems,  wave  number  plays  a. 1  .he- 

=iSSSlH=SH“;=s 
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Koitcr's  and  Carrier’s  justi.ication  of  the  use  of  substitute  kernels  in  the  transform 
plane  has  the  immediate  corolLty  of  validating  the  use  of  simplified  boundary 
conditions  to  describe  the  physics,  since  the  kernel  is  in  effect  the  transform  realiza¬ 
tion  of  the  boundary  condition  operators.  Our  use  of  the  set  (12)  as  boundary 
conditions  is  however  a  preferable  procedure  for  it  is  more  direct  and  we  can  justify 
it  on  physical  grounds  by  using  the  Haskell  method  backwards.  That  is,  we  can 
transform  the  conditions  in  (12)  from  the  crust  -mantle  interface  to  the  free  surface. 


tw 

FlO.  4.  (a)  The  crustal  layering  on  either  side  of  the  land  sea  interface  can  be 
thought  of  as  a  continuous  distribution  of  shunt  transmission  line  networks 
on  cither  side  of  the  discontinuity  if  the  exact  Haskcll-Thomson  matrix 
formulation  is  used,  (b)  The  simplified  procedure  is  analogous  to  replacing 
the  transmission  line  shunts  by  lumped  parameter  equivalents  which  have 
similar  characteristics  in  the  relevant  region  of  the  wave  number  plane. 


Naturally,  they  will  not  imply  that  the  surface  stresses  vanish  except  for  Rayleigh 
waves.  However,  any  incident  ray  whose  wave  number  is  close  to  the  Rayleigh  wave 
number  will  not  induce  'ippreciable  stresses  at  the  surface.  Indeed,  for  P  or  Sf' 
waves  striking  the  free  surface  at  grazing  angles  of  incidence,  the  equations  in  (12) 
describe  the  physics  exactly.  Essentially  all  the  diffracted  rays  arc  tangential  to 
the  surface  discontinuity  and  for  waves  of  this  type  the  stress  mismatch  at  the  surface 
will  be  small.  In  any  event,  more  terms  can  be  added  to  the  simplified  boundary 
conditions  after  the  fashion  of  the  analysis  in  (7,8,9)  to  get  any  desired  degree  of 
surface  matching. 

4.  Solution 

In  terms  of  *p(x,y)  !’.nd  Mx,  y)  the  synthetic  boundary  conditions  (12)  at  the 
crust-mantle  interface  can  be  written 

o.  jt-o-  (14) 
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is<p+f. “0*  y  “  o 

1  *.t 

(15) 

^0  km  i  ,  A 

rx+r~*‘-  “•  '-0 

(16) 

S  =  +kK  sin  0, 

(17) 

where 


is  a  coefficient  that  always  involves  kR  owing  to  the  suppressed  z-variation. 

Let  R(«)  denote  the  hermitian  matrix  (rtJ)  where  the  coefficients  r(J  are  identified 
as  polynomial  components  of  the  realization  of  the  differential  operations  (14—17)  in 
the  transform  plane  when  x  is  positive,  viz. 

/** r„  -IS  -la  \ 

«K*)-(  IS  fca/r*  0  (18) 

\  i«  0  kg/TnJ 

Likewise,  let  L(a)  be  the  matrix  counterpart  of  (19)  but  with  parameters  kL  and  Ft 
appropriate  for  the  left  boundary.  The  right  matrix  R  is  singular  when 

a  «*  ±r  =  ±k*cjs©,  (19) 

;nd  the  left  matrix  L  is  singular  when 

a  ~  ±1  =  ±(kl~kii in2©)*.  (20) 

These  points  correspond  to  the  geometric  acoustics  poles  for  the  right  and  I  *) 
boundaries  respectively.  Tne  relationships  (19)  and  (20)  are  seen  to  define  the 
correct  angles  of  refraction  at  the  interface. 

We  introduce  the  unknown  vector 


and  a  constant  vector  m  describing  the  mismatched  components  of  the  incident  wave 
upon  the  left  boundary 


where 

kKr h 

"sMWIV^sine  (23) 

fMj  =  KMVTt-Mcose 

A.  Function-theoretic  method.  Vajnshtcjn  (13).  Karp  (14),  and  Clemmow  (15)  have 
introduced  a  rather  general  method  for  solving  dual  integral  equations.  The  VKC 
procedure  can  be  generalized  to  solve  vector  problems  with  either  the  exact  or 
simplified  boundary  conditions. 
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A  plus  ( or  minus)  superscript  on  a  matrix  or  vector  quantity  will  indicate  that 
all  of  its  components  arc  analytic  in  the  upper  (or  lower)  a-pianc  above  (or  below) 
the  contour  of  integration  <€.  Unless  otherwise  specified,  tnesc  quantities  vanish 
uniformly  at  infinity  in  their  rclevan'  haT  plane  of  regularity,  By  Jordan  s  lemma, 
the  boundary  condition  on  the  right  wifi  be  satisfied  provided  that 

R(a)x  -=  f  *  (at)  t24) 

where  i;»  ir  unknown  ar  .A  On  the  left  boundary  we  need  satisfy 

Ua)x  f—  «  p'(«).  @5- 

*4  t 

The  term  mM-rk.)  represents  a  collection  of  poles  whose  residues  contributions 
represent  the  driving  e'ect  of  the  incident  Rayleigh  wave  upon  the  left  boundary. 
In  order  to  solve  these  equations,  we  substitute  (24)  into  (25)  and  obtain 


K(e)f(«)+; 


8'(a) 


where  the  matrix  to  be  ‘split’  is 

K(«)  -  MaXRWr1.  <27> 

For  the  rimplified  boundary  conditions  ti.e  entries  in  the  matrix  kernel  IC(a)  art 
t0  be  rational  functions  of  the  transform  variable  a  so  that  bv  elementary 

means  we  can  split  K(a)  as 

K(a)  =  (N"(a))“‘P*(«)  <28> 

For  the  exact  formulation,  however,  contour  integration  would  be '  *® 

perl  >rm  the  factorization.  In  any  event,  the  matrix  factors  P  («)  and  M  W  are 
regu  ar  in  their  respective  half  planes,  however  unlike  the  vectors  f  (a)  or  0  (u/ 
they  do  not  vanish  at  infinity  but  approach  a  constant  value. 

If  we  add  and  subtract  the  same  term,  we  can  rearrange  equation  (26)  into  the 


P*(«)fV)4- 


W  (r^  =  N'(x)3'(a)+-r-lN'(-r)-hl‘(*)]m- 
*4-  r  *  +  f 


Each  side  of  this  equation  i,  regular  in  a  half  plane  and 

tion  of  an  nuire  function.  Owing  to  the  assumed  growth  off  W  « »nd ,  a l  («)  this 
entire  fun  ion  vanishes  uniformly  at  infinity  and  must  thus  be  the  null  constant 
Hence  we  can  set  each  side  of  (29)  to  zero  and  obtain 

.IP*(a)]-,N-(-r)m 

W  a  +  r  (30) 

«-(<)«  _EM-!-tN-(-r)-N-(a)lm. 

*  '  «4-r 

With  the  aid  of  the  first  of  these  expressions,  the  desired  solution  x  is  seen  tc  be 

IRfair^WPW-t-rlm  '  (31) 


a4-r 
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B.  Algebraic  method.  For  rational  kernels  the  general  procedure  described  in 
section  A  is  rather  abstruse  and  a  more  direct  approach  can  yield  the  solution  in  a 
simple  algebraic  fashion;  at  the  same  time  the  transparency  of  the  analysis  will 
offer  an  insight  into  the  mathematics  that  will  suggest  a  generalization  that  can  be 
used  to  solve  more  complex  problems. 

The  matrices  L  and  R  are  sue*  that  the  only  singularities  we  can  have  in  the 
solution  x  are  poles  in  the  complex  plane  at  a  «  ±r,±l.  We  know  that  the  vector 
x  must  have  entries  which  are  rational  functions  of  a.  The  poles  of  x  can  only 
arise  at  a  «  ±r  (incident  and  reflected  waves),  and  ?*.««■  —I  (transmitted  wave) 
for  the  residue  wave  at  a  =  +  /  is  a  non-physical  contribution.  Hence,  x  must  be 
inversely  proportional  to  (a  —  r)  (a  +  r)  (a  +  0.  ar*d  we  can  immediately  concentrate 
attention  on  the  numerators  of  its  components.  That  is,  x  must  be  of  dte  form 

fl,  +  fli(a-r) 

sin  0(a4  +  <ij(a— r)(a 

<r„  cos0[a,  +  fl}(a- 


1 


(a-r)(a+r)(a+0 


iT 


where  the  a,  arc  five  constants  to  be  determined.  The  order  of  the  polynomial 
..-.tries  in  x  is  easily  established  if  we  anticipate  that  we  need  use  Jordan’s  lemma 
after  multiplication  by  either  the  operator  R  or  L  The  five  unknown  constants 
a,  can  and  must  be  so  chosen  that  all  boundary  conditions  are  obeyed.  For  x 
positive,  we  need  satisfy 


-L  fRxe^dor-O. 
2ni  J 

* 


x  >  0 


(33) 


and  the  form  of  (32)  is  such  that  it  automatically  meets  this  constraint,  since  the 
only  residue  contribution  in  the  upper  half  plane  is  a  reflected  Rayleigh  wave. 

For  the  boundary  condition  on  the  left  side  we  need  satisfy 


bit; 


Lx+~'Vl"da  ■*  0,  x<0 


a+rj 


(34) 


and  by  Jordan’s  lemma  this  requirement  will  be  met  if  the  residue  coefficients  at 
a  «  -r,  -I  are  null  That  is  wed  satisfy 

(«+r)L(-r)x(-r)+m  =  0,  (35) 


and 


(«+/)L(-f)3(-/)  =  0. 


(36) 


Equation  (36)  is  a  homogeneous  relation  and  can  be  satisfied  by  specifying  two 
constants  to  make  th*  residue  contribution  proportional  to  a  transmitted  Rayleigh 
wave  that  matches  the  left-hand  conditions,  viz., 


T,  sin  ©[a,  +  a,(l+ r)(f + *=  —  rt sin  ©.(a,- <j,(/+ r)) 

-rj,cos0[<J|-flj(/+r)]  =  +rtcos©l[flI-fli(/+r)). 


The  three  equations  in  (35)  represent  the  left  boundary's  reaction  »hat  cancels  the 
incident  fields  4nd  together  with  th<  two  in  (37)  wc  can  determine  the  five  constants 


f 
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a,  uniquely.  A  simple  transformation  will  have  the  result  of  partitioning  these 
equations  into  easily  solvable  forms.  Let  i/j  **  T ijOj  where  the  matrix  Ty  is 

(rt  -2rrt  o  o  o  \ 

r,  sin  6  0  2r,rJ,sin©  2rr*sir,©  0  \ 

r*cos©  0  0  0  -2rr, cos 0  I  (38) 

rt  -(/+r)  rt  o  o  o  I 

i  r,  o  Kf+or,  (f+r)r„  r  / 

In  terns  of  the  bjs,  the  equations  (35)  and  (36)  can  be  expressed  by  the  matrix 
relation 


-S  -r 


-S  kL 


0  sin  0,  —  sin  ©/  \  6« 
Pi  Pi  /  \by 


(c2  sin  0 
cx  cos© 


where  the  parameters  p„  p2,  and  p3  are 

I+r  c. 

p, -  — T-r*cose. 

A 

p,  -  -£r..  W 

pj  *  rtcose,+rnco$©p 

The  first  three  unknowns  6„  and  6,  are  easily  found  to  be 

h,  «  rL, 

b2  =  ir(/-r)r„sin©,  (40 

ft,  =  r,cos0 

and  »hen  from  the  last  two  equations  in  the  matrix  relation  (39)  we  can  calculate 
the  remaining  parameters 

ziy^cosoq’-r*) 

4 "  r«cos©+rtcos©,  (42) 

,  sin©,. 


The  inverse  transformation,  a}  =  Ty  xbt,  determines  the  ajj.  In  particular,  we  shall 
require  expressions  for  at  and  a2 


...  2r,cos©  1 

a,  «  2r(l  r)[l-rj(COS0  +  rtCOSe,J 

2(l+r)rK  cos© 
rj,cos8+rtcos©, 


O]  -  2 r— 


since  these  numbers  arc  required  to  evaluate  the  residue  cocfTicicnts  —  at/I2r(/  +  r)J, 
and  [a,— a2(/+ —  ra J  at  the  poles  a  »=  r  and  —  /  respectively. 

5.  Discussion  of  the  solution 

The  elementary  character  of  the  solution  (31)  or  (36)  permits  a  ready  evaluation 
of  the  scattered  fields.  For  y  **  0  and  x  positive  the  entire  field  is  a  collection  of 
residue  contributions  which  can  be  grouped  into  an  incident  Rayleigh  wave,  and 
a  reflected  Rayleigh  wave  with  a  relative  ..nplitude 

a  »  C3S0 Z  (r cA  1  -(P| Jv*)[  Q)  r44v 

cos  0  -f (I'j/r *),/[ I  -  (oJvK)1  sin*  ©]• 

Our  con  "ntion  will  be  to  describe  the  coefficient  of  the  comprcssional  potential 
of  the  Rayleigh  wave  as  its  amplitude.  Although  our  solution  includes  scattered 
body  waves,  their  contribution  is  null  on  the  reference  plane  owing  to  our 
formulation. 

On  the  left  side  of  the  boundary,  the  incident  Rayleigh  wave  is  cancelled  by 
pole  contributions  as  it  should.  The  remaining  residue  terms  are  seen  to  represent 


■too  -so  o  so  too 
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Fie.  3.  Contours  of  constant  values  of  the  transmission  coefficient  f  for  an 
incident  Rayleigh  wave  at  normal  incidence  f0  •*  0)  upon  crustal  dis¬ 
continuity.  For  this  special  case  there  is  no  variation  of  j  with  respect  to 
the  velocity  contrast 


V«(«  witt  SifTrullo*  it  crwstit  totawUwwttks,  t 


a  transmitted  Rayleigh  wave  of  amplitude 

_  2  cos  0 

+  “  cos©+(rt/r«K/H  -{vjv^  sin*  0]' 


This  wave  travels  in  a  direction  0,  which  is  given  by 

_  vt  ■  ~ 

sin  0,  «  —  sin  0, 


(45) 

(46) 


a  familiar  expression  for  the  refraction  of  a  ray  at  an  Interfax.  In  general  both 
phase  velocities  vL  and  p*  depend  upon  frequency,  and  as  a  consequence,  the  angle 
0,  is  frequency  dependent.  Thus  the  diffraction  of  a  Rayleigh  pulse  incident 
obliquely  upon  a  crustal  interface  will  sort  the  transmitted  spectral  components 
into  different  directions  after  the  fashion  of  an  optical  prism.  In  addition,  if  0 
exceeds  a  critical  angle 

0, «  sin“,(pt/pt)  (47) 


then  total  reflection  takes  place.  Just  as  in  optics  or  acoustics,  this  can  only  happen 
if  the  wave  is  incident  from  a  slow  to  a  fast  medium.  Finally,  we  note  that  if. 
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Fm.  6.  Contours  of  constant  values  of  the  transmission  coefficient  f  for  an 
*  incident  Rayleigh  wave  at  oblique  incidence  (9  -  30*)  upon  a  crustal 

discontinuity. 


-wo 


so 


50  0 
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Fig.  7.  Contours  of  constant  values  of  the  transmission  coefficient  S'  for  an 
incident  Rayleigh  wave  at  oblique  incidence  (0  »  60*)  upon  a  crustal 
discontinuity.  Note  that  total  reflection  takes  place  if  cw'vL  exceeds  2/,/J. 


vjv„  *  1,  then  there  is  no  angular  variation  in  31  and  S'  for  a  discontinuity  in  T. 
Figs.  5-7  illustrate  the  variation  of  the  transmission  coefficient  as  a  function  of  the 
relative  discontinuities  in  the  Rayleigh  wave  parameters  for  fixed  angle  of  incidence. 
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TEl.KSKJSMIC  RESPONSE  Of  A  UNIFORM  DIPPING  CRUST 

(PART  I  OF  A  SERIFS  ON  CRUSTAL  EQUALIZATION  OF 
SEISMIC  ARRAYS) 

By  Julius  Kane* 


ABSTUACf 

Mo*t  seijm'u.  theory  it  confined  to  o  considerotion  of  cruitol  ttructurei  thot  con  be 
obttrocted  01  plone  poroltel  eloitic  loyer;.  For  such  configurotion*,  the  reipome 
of  eoch  element  of  o  leismic  orroy  will  be  timilor  to  ony  other  except  for  o  time 
deloy.  Signol  enhoncement  of  on  orroy  located  on  tuch  an  idealized  structure 

con  be  accomplished  by  suitable  time  deloyt  of  the  individual  tracei  followed 

* 

by  superposition. 

On  the  other  hond,  if  the  erustol  structure  is  onything  but  a  plone  parallel 
configuration,  the  signols  received  by  the  individual  elements  will  not  be  identical 
to  one  another,  but  will  include  distortion  effects  chorocterisFc  of  the  tocol 
geometry.  As  o  result,  the  records  of  seismic  arrays  located  on  realistic  crustal 
configurations  will  hove  to  be  equalized  to  some  standard  reference  if  optimum 
signal  processing  is  to  be  achieved. 

In  this  paper  w*  introduce  a  ray  procedure  for  the  calculation  of  theoretical 
seismogroms  for  the  teleseismic  response  of  an  array  of  stations  located  above 
a  uniform  dipping  crust  (wedge-shaped).  In  terms  of  this  mathematical  model,  we 
demonstrate  the  signol  distortion  effects  of  the  geometry  and  discuss  equaliza¬ 
tion  techniques  that  will  permit  a  superior  recovery  of  the  desired  signal. 

iNTHOUUCTtON 

The  outcome  of  any  experiment  is  subject  to  random  errors,  and  otic  way  of 
reducing  uncertainty  Is  to  rc(>cat  the  measurement  many  times.  Before  the  results 
of  the  individual  measurements  are  averaged  there  is  one  canal:  one  must  be  sure 
that  the  same  experiment  has  been  performed-  In  seismology,  the  need  for  enhancing 
low  level  signals  relative  to  background  noise  is  aerte,  and  for  this  reason  elaborate 
arrays  are  being  constructed  to  monitor  telcscisms.  Owing  to  various  reasons  it  is 
not  possible  to  position  the  site  of  each  station  or  element  ill  the  array  on  terrain 
identical  to  that  of  its  neighbors.  As  a  result,  each  element  will  respond  to  the  same 
signal  in  a  perturbed  fashion  characteristic  of  the  local  crustal  structure.  This 
principle  is  of  course  well  known,  and  a  study  of  the  signal  variations  across  an 
array  of  closely  spaced  stations  for  ground  shots  is  one  way  of  mapping  subterranean 
strata.  In  source  mechanism  studio  of  sciatic  sources,  the  situation  is  somewhat 
different  from  the  field  mapping  environment.  The  vibrations  can  be  considered  to 
be  that  of  a  plane  wave,  the  signal  arrives  from  below  rather  than  above,  and  the 
spacing  of  stations  is  comparatively  wide  as  compared  to  those  used  in  field  surveys. 
Nonetheless,  for  a  varying  crustal  structure,  the  response  of  any  element  will  be 
distorted  by  terrain  effects  so  that  each  station's  signal  should  be  equalized  to  some 

*  Present  address:  Space  Science  Center,  University  of  California,  Los  Angeles,  California. 
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standard  reference  Iwfnrc  processing  to  reduce  noise  effects.  Severn  I  questions  arc 
pertinent:  how  severe  might  this  effect  lx?  for  real  in  tie  crustal  configurations,  and  If 
the  effect  is  important,  hmv  ran  the  signals  h'  processed  so  us  to  remove  the  effects 
of  local  geometry?  In  this  paper  we  consider  this,  first  question  by  determining  the 


Fro.  I.  Au  incident  plane  wave  arrives  at  a  positive  angle  0  measurer!  with  res|»cct  to  the 
normal  to  the  free  surface,  and  illuminates  a  crustal  wedge.  An  array  of  stations  A,  R,  C,  are 
located  on  the  free  surface;  their  images  are  shown  at  A’,  .1*.  A",  etc.  The  ray  paths  for  the 
velocity  variation  given  in  Tnide  l  are  shown;  the  T-patlis  are  drawn  ns  full  lines,  and  the  S- 
paths  as  dashed  lines  Any  ray  which  will  fail  to  intercept  another  im  ,gc  plane  is  riot  drawn. 

theoretical  seismograms  that  an  array  of  stations  might  receive  if  located  on  a  dip¬ 
ping  crustal  structure.  The  advantage  of  this  particular  choice  of  geometry  Is  that  it 
allows  a  rigorous  calculation  of  the  dominant  portion  of  the  scismogrrtns  for  any 
assumed  variation  of  the  incident  wave. 

In  addition  to  the  theoretical  significance  of  the  present  study,  the  analysis  has 
practical  importance  for  preliminary  analyst*  at  the  Caltech  Sctsmological  Lab- 


TF.I.ESEInMIC  HESPONSK  or  VMKOHU  OUTING  ClIVhT 


oratory  (Xuui,  196,*>)  mid  the  U.S.  Geological  Survey  (Roller  ft  at,  1904)  indicates 
that  the  Touto  Forest  array  is  located  on  a  dipping  crustal  structure  with  a  dip  of 
approximately  S°  and  an  estimated  strike  of  N  17°W. 


Formulation 

We  consider  au  array  of  stations  indicated  by  A,  B,  and  C  located  on  a  dipping 
crustal  structure  (Figure  1).  An  incident  P-  or  5- wave  arrives  at  an  arbitrary  angle 
of  incidence  and  illuminates  the  configuration.  Our  task  is  to  determine  the  subse¬ 
quent  ground  motion  for  a  given  time  variation  in  the  inconvng  wave. 

The  task  posed  is  an  intractable  boundary  value  problem  a  the  entire  solution  is 
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to  Ire  found.  However,  tin*  fu-lil  rim  Ik*  separated  into  tlto  sum  of  two  components:  a 
multiply  reflected  field,  mid  :m  involved  cylindrical!}'  sputtered  fielil  rc-rudiating 
from  n  virtual  source  at  the  vertex.  Fortunately,  it  run  he  shown  that  this  compli¬ 
cated  ihlfraetion  held  does  tint  eoutriliiite  significantly  to  the  ground  motiott  for 
several  reasons.  First,  the  major  lobe  of  this  rylimlrirul  field  is  direeteil  downward 
away  from  the  stations  where  it  smooth*  the  almtpl  transition  between  the  different 
reflected  fields  from  either  boundary.  Second,  it  will  tie  shown  that  the  vertex 
intercepts  only  an  infinitesimal  fraction  ol  the  rays  that  enter  the  crustal  wedge. 
Indeed,  a  realistic  crustal  configuration  may  nol  have  a  vertex  since  the  actual 
crustal  structure  may  lie  wedge  sha|K*d  only  in  the  immediate  vicin.ty  of  the  array. 
Third,  the  cylindrical  field  decays  with  distance  s«.i  its  amplitude  falls  off  rapidly, 
mid  its  contribution  at  the  stations  will  be  negligible  since  we  shall  assume  that  they 
are  relatively  distant  from  the  discontinuity.  For  these  reasons  wc  shall  neglect  this 
cylindrical  field,  foi  while  its  determination  would  be  an  interesting  (and  difficult) 
mathematical  task,  this  effort  would  be  ol  limited  geophysical  interest. 

In  the  time  domain,  the  seismograms  can  be  determined  by  following  the  paths 
of  the  multiply  reflected  anil  refracted  rays  within  the  wedge.  Unlike  the  case  of  a 
plane  parallel  layer,  the  ray  tracing  procedure  can  be  very  tedious  unless  it  is  done 
properly.  A  very  convenient  procedure  is  to  reflect  the  wedge  rather  than  the  rays 
as  indicated  in  Figure  1.  In  this  diagram  the  stations  arc  multiply  imaged  at  A\  B , 
(f,  •  •  •  .4",  B",  C",  •  ••  .4”,  /?".  O",  ■  ■  ■  etc.  In  this  diagram  the  path  of  a  multiply 
reflected  P-  or  6'- wave  is  simply  urawu  as  a  straight  line,  since  the  angle  of  incidence 
equals  the  angle  of  reflection.  However,  the  paths  of  rays  following  P/S  or  S/P 
conversions  are  complex  broken  lines.  To  compute  the  seismogram  wc  consider 
each  station  ami  all  its  images  as  one  and  combine  the  outputs  of  each  imaged 
station  on  one  seismogram  following  the  laws  of  vector  addition. 

The  lice  diagram  hns  several  advantages.  For  otic  thing,  it  clearly  shows  that 
any  ray  that  enters  the  crustal  wedge  other  than  at  the  origin  trill  avoid  the  vertex  and 
eventually  migrate  to  infinity.  Tims  the  edge  intercepts  only  an  infinitesimal  fraction 
of  the  energy.  Another  advantage  of  tin  presentation  is  that  it  shows  that  the  path 
of  any  my  has  a  counterclockwise  drift  and  it  w  ill  cycle  from  one  layer  to  the  next 
in  a  rc|H“tilivc  fashion.  This  means  that  all  the  calculations  can  easily  be  done  in 
the  computer  by  reiterating  a  two-stage  computer  loop,  which  describes  the  two 
reflecting  surfaces.  An  additional  advantage  of  the  tree  diagram  is  that  it  permits 
visual  interpret  at  ion  of  a  vast  amount  of  computer  output  which  gicutlv  assists 
the  computer  debugging  0|K*r:itiott.  In  Figure  ?  wc  give  examples  of  computcr- 
drawii  tree  diagrams  for  a  selection  of  various  angles  of  incidence  U|ion  a  crustal 
wedge. 

lit  all  our  illustrations,  wc  have  chosen  the  parameters  in  T'able  1  for  calculation 
of  the  reflection  coefficients.  The  values  used  arc  proportional  to  the  values  shown. 

With  pm|)er  interpretation  the  program  used  to  generate  the  tree  diagram  contains 
sufficient  information  to  const rurl  a  theoretical  seismogram  for  a  station  arbitrarily 

1  Our  data  is  ended  on  magnetic  tape  by  an  HIM  7000.  and  then  irnnscriticiS  into  visual  out- 
pui  by  the  Siromberg  Carlson  tO’JO  which’  nulomaltcnllv  drives  and  photographs  the  face  of  a 
special  Chartclroii'*  oscilloseo|>c  rtil>c. 
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TAUI.K  1 

AOOPIKD  VELOCITIES  AND  DENSITIES 
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located  on  the  wedge  surface.  The  diagram  is  diav.ii  hy  the  computer  with  the 
assumption  that  the  incident  ray  strikes  the  lower  wedge  urfaic  at  any  (mint  say 
unit  distance  from  the  vertex.  This  jxhnt  locates  the  trunk  of  the  tree,  and  the 
brunches  spread  out  according  to  the  laws  of  reflection  and  refraction.  Each  branch 
that  intersects  one  of  the  image  planes  carries  with  it  nine  pieces  of  datu  whi  n 
combine  to  form  a  vector  E(c, ,  Cj ,  •  •  •  n)  called  an  emit.  Each  event  represents  an 
incident  plane  wave  striking  one  of  the  wedge  surfaces  tog?tlicr  with  two  reflected 
rays.  The  components  of  this  nine-dimensional  vector  are  identified  ns  follows: 

(1)  Bounce  number.  The  first  coordinate  cj  will  be  referred  to  as  the  liouncc  num¬ 
ber.  The  magnitude  of  the  bounce  number  fi  locates  the  image  plane,  and  its  parity 
designates  whether  the  event  corresponds  to  an  incident  shear  (negative),  or  com- 
prcssionnl  (positive)  ray.  The  bounce  number  has  any  value  dfcjn  where  n  is  any 
integer  other  than  zero.  At  the  low  cst  surface,  tire  crust -mantle  interfare,  the  bounce 
number  Ci  has  the  value  +J  if  the  excitation  is  a  comprcssional  wave.  An  integral 
bounce  number  indicates  that  the  ray  is  incident  upon  the  free  surface  and  this 
associated  event  will  contribute  V  the  theoretical  seismogram;  otherwise  the  ray 
is  incident  ujiou  the  clastic  interface  and  will  not  contribute  to  the  resjiousc  until 
another  reflection  takes  place 

(f )  Travel  lime.  The  second  coordinate  is  the  travel  time  required  for  the  ray 
to  reaeh  the  relevant  image  plane  from  the  root  of  the  tree.  At  the  base  of  the  diagram 
we  set  Cj  *»  0,  and  the  travel  time  will  accumulate  anthmet.cally  from  this  origin 
as  the  ray  progress  from  bounce  level  to  bounce  level. 

( 3 )  Radial  distance.  The  third  coordinate  is  the  radial  distance  from  the  wedge 
vertex  to  the  ray  intersection  at  each  bounce  level. 

(i)  Incident  magnitude.  The  fourth  coordinate  c»  is  the  relative  magnitude  cf 
the  incident  ray  as  it  strikes  the  image  plane  indexed  by  r, .  At  the  root  of  the  tree 
we  shall  normalize  <•«  =  1,  and  this  coordinate  will  progress  geometrically  from  level 
to  level  in  rcsjiousc  to  a  factor  proportional  to  the  appropriate  reflection  cocflicient. 
For  convenience,  and  greater  realism,  « lie  reflection  coefficient  used  at  the  earth’s 
surfurc  was  not  that  appropriate  to  the  free  surface  of  an  elastic  solid,  but  rather 
that  for  a  solid-fluid  interface  wherein  the  elastic  parameters  for  air  arc  used. 

(5)  Sine  of  (he  angle  of  incidence  The  fifth  coordinate  c4  =  sm  d  where  d  is  the 
angle  of  incidence  of  the  incident  rav.  It  is  preferable  to  use  the  sine  of  d  rather  than 
d  itself  as  a  variable  for  the  »aws  of  reflect  icn  and  refraction  depend  algebraically  ujxm 
-  sin  d  and  in  a  transcendental  fashion  upon  d.  By  use  of  this  variable,  we  ob- 
tain  the  advantage  of  mnrkedly  accelerating  the  computer  program  for  the  time- 
consuming  trigonometric  subroutines  need  never  be  used  in  the  iterative  loops. 
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""  mr“"ml  l"'l«  *»  «»  vcr.i,  ,1  acrtrf  „t  the  eurtnit 

(«)  llrjlrclol  I'.uwinilmle.  The  -'.Ml.  coonlinalc  e.  i,  the  magnitude  of  .he  re- 
fleetal  roiii|ire.uno:i:il  my.  I„  eah-f.latieg  «  haVC  „„ ,  o( 

coeflie'cuts  ,  e|«.,,,h„R  the  |«ri,r  of  the  h.eidc.t  my,  lvhclhor  ,hc  Mice! 

rS':;1'  •*  or  “  *■* 


RI,P( c5),  /?/,,(r»), 

wherein  «hc  second  letter  indicates  whether  the  reflection  takes  place  at  solid-flu, d 
imerfarc  (/?F),  or  at  a  solid-solid  interface  (/?/),  the  subscript  8p  for  example  in- 

dieatcs  that  the  reflet,  on  is  from  an  incident  shear  wave  In  a  reflected  comprcioual 
wave.  1 

(7)  Sine  of  the  angle  of  the  reflected  P-umc  The  coordinate  e7  is  the  sine  of  the 
reflected  P- wave  as  computed  by  Snell’s  law, 

and  has  a  sign  opposite  to  the  op|>ositc  sign  to  the  incident  w..ce 
(«)  Reflected  S-magnitudc.  The  eighth  -oordiuatc  e.  is  the  magnitude  of  the 
reflected  shear  ray  and  requires  the  proper  use  of  one  of  four  expressions  for  its 
calculation;  these  arc  written  in  an  obvious  notation  as 


**«<*>,  Rf'p>(et), 

(9)  Sine  of  the  angle  of  the  reflected  S-iratc.  The  last  coordinate  c,  is  the  sine  of 
the  reflected  5-wavc  as  computed  by  Snell’s  law,  and  has  the  opposite  sign  to  the 
incident  wave.  b 

ITic calculation  of  the  tree  diagram  pr.ice.sls  in  a  very  simple  and  straight -forward 
manner  once  the  geometry,  elastic  parameters,  and  excitation  are  specified.  Typi¬ 
cally,  the  cxeita’ion  will  be  chosen  as  a  P- wave  incident  from  below  the  wedge,  and 
the  time  variation  of  the  stress  along  the  incident  ray  will  be  assumed  to  be  given 
Let  A  be  the  relative  amplitude  of  this  ray.  l>u  striking  the  elastic  u,  Jge  two 
rays  w,U  pc  >ctratc  the  crtstnl  structure,  and  we  have  the  following  initial  vector 
“  '*»  >  *• »  ‘  * '  f» )  to  start  the  iterat  ion  scheme  wherein 

<«*  =  +i, 

e *  —  1. 

e%  -  sin  (0  -  a), 

c*  =«  -p;  sin  (©  -  o), 

l> 

t*  ■  ^7  sin  (0  —  a). 

«V 

in  these  c.  pressious,  0  is  the  angle  of  incidcu.  e  of  the  wave  measured  with  respect 
to  a  vertical  on  the  free  surface  (Figure  1),  and  is  th  vedge  angle,  and  r/  and  r/ 


ct  =  0, 
c*  *=  *4, 

te  -  dr„[sin  (0  -  «)],  (1) 

c8  =  .ITp.Km  (0  -  a)]t 
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are  the  velocities  of  the  shear  and  roiupicsMomd  waves  in  the  uutnHc.  Within  the 
crustal  wedge,  the  cxeitntioa  is  split  into  two  rays  with  particle  amplitudes  projtor- 
tiounl  to  T„  and  T ,, ,  the  transmission  eoeflieients  for  a  plane  wnvc  incident  upon 
nn  clnstic  interface.  This  initial  vector  is  used  to  stati  the  iteration  scheme.  The 
logic  behind  the  eompntntional  process  is  cosily  described  in  two  steps. 

Step  1.  At  nny  image  level  suppose  all  the  events  have  been  calculated,  livery 
event  is  placed  in  temporary  storage,  and  in  addition  if  ct  is  an  integer,  then  the 
corresponding  event  is  placed  either  in  permanent  storage  or  written  ott  tape.  In 
partirulnr  (1)  would  be  the  only  event  for  the  lowest  level.  In  general,  if  «i  ■  Jr, 
then  there  will  be  a  maximum  of  2'"1  events  to  be  calculated  at  that  level.  In 
practice,  this  number  will  1m?  substantially  less  for  we  shall  discard  any  incident 
ray  unless  its  relative  magnitude  exceeds  some  threshold  value. 


Fio  3.  The  geometry  behind  the  deration  scheme. 


Step  2.  Once  all  the  events  appropriate  to  an  image  level  have  been  calculated 
and  stored,  then  we  use  each  event  in  temporary  storage  to  calculate  a  new  pair  of 
events  for  the  next  image  plane.  The  calculations  at  this  stage  arc  essentially  repeti¬ 
tions  of  one  basic  series  of  algebraic  operations.  I-et  E\e i,  ef,  •  •  •  r,‘)  be  a  typical 
event  in  temporary  storage.  The  two  new  vectors  to  be  generated  can  be  denoted 
asE*  and  £* since  they  correspond  to  the  reflected  P-  andSF-waves.  T  he  components 
of  'hese  new  vectors  will  be  algebraic  functions  of  the  preceding  coordinates. 

Let  us  introduce  the  abbreviations  for  the  constants 

C  *  cos  a,  S  »  sin  a,  (2) 

.  and  designate  ihe  radial  distances  r,  and  r,  shown  in  Figure  3  simply  as  r 
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fashion  until  the  data  is  depleted.  Once  the  event*  in  temporary  storage  arc  dc- 
plotcd  we  cycle  back  to  stage  1. 

The  Tiikoheticai,  Seihiiogham 

A.  Preliminaries.  In  the  fashion  just  described,  we  can  accumulate  a  rise  voir  of 
events  for  any  desired  configuration  of  system  parameters.  However,  these  vectors 
have  little  value  until  they  are  converted  into  a  theoretical  seismogram.  In  the 
conversion  process  a  great  deal  of  simplification  is  introduced  if  we  remind  ourselves 
that  the  geometry  of  the  wedge  problem  is  such  that  U  is  specified  completely  by  angles, 
and  that  it  incorporates  no  fundamental  unit  of  length.  This  fortunate  circumstance 


Fio.  4.  The  scaling  of  any  branch  on  the  tree  diagram  proceeds  by  similar  triancles 
and  a  causal  correction  g-ves  the-  arrival  time  of  any  event  at  any  station 

will  permit  us  to  use  the  scaling  r  i  similar  triangles  to  normalize  each  event  to  the 
location  of  any  station. 

^rc  havc  calculated  the  tree  diagram  with  its  base  at  r  =  1,  and  let  the  branches 
of  the  tree  fall  wbarc  they  may.  Suppose  we  have  a  station  A  at  r  ■»  a,  then  we 
need  to  know  the  travel  times  when  the  events  reach  r  -  a  at  the  free  surface  or 
one  of  its  images.  In  general,  we  need  alter  the  travel  time  e,  and  radial  distance  e, 
so  that  the  event  is  properly  normalized  to  give  the  correct  contribution  at  any 
station.  By  similar  triangles  (Figure  4),  we  sec  that 

•r  _  ?’e*  _  0  “  a/e,)  sin  (9  —  a) 

e,  p/  .  (*) 

is  the  corrected  travel  time  for  the  initial  occurrence  of  the  event  at  station  A. 
The  second  term  in  this  expression  represents  a  causal  correction  for  the  excitation 
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wave  arrives  at  different  parts  of  (hr  ba-o  d  (he  crustal  wedge  at  different  times. 
\\  e  note  that  the  corrected  travel  time  r„  depends  upon  the  ratio  a/e j  :  us u  rule,  this 
fraction  increases  without  limit,  with  increasing  bounce  number  or  higher  older 
revertwiiit  ions.  As  a  result ,  tin  nj  a  station's  res  j, rase  relative  to  another  at  a 

dUfercnt  radial  distance  ineiriisr  udh.  tune.  In  other  words,  one  effect  of  a  dipping 
crustal  structure  can  l>e  given  the  interpretation  that  cadi  component  of  an  army  o* 
stations  will  res|>ond  to  the  same  signal  as  though  it  were  written  on  a  rubber  sheet. 
The  stretching  of  the  s.teet  will  be  piojKn  t  i<*iml  to  both  time  and  the  relative  dirtnucc 
of  the  station  from  the  wedge  vertex. 

B.  The  Ground  M otion  IJicory.  If  wc  neglect  pure  distortions,  for  two  dimensional 
motion,  the  displacement  7  «  (u,  v)  of  an  elastic  solid  is  given  by 


difi 

u  “  r-  +  r ~ , 

dx  dy 


dy  dx ' 

where?  and  are  solutions  of  the  wave  cijuations 


vV  -  - 

‘V* 


1  dV> 


dr- 


0, 


(7) 


(8) 


In  our  analysis,  wc  have  confined  our  attention  to  the  ray  portion  of  the  seismic 
field  so  that  ?  and  ^  will  be  plane  waves  of  the  pat  t  irular  form 


V  -  /  ( t  ±  —  ."in  i?,  ±  —  cos  d. ) 

\  l’p  vr  } 

~  /  (^l  ±~  sill  !?,  rb  ^  cos 


(9) 


where /(f)  is  an  arbitrary  function,  of  t  and  «?p  ,  it,  are  any  parameters.  An  important 
feature  of  plane  wave  reflection  and  refraction  is  that  the  functional  form  of  /(/)  is 
common  to  both  incident  and  reflected  waves  r.  long  as  the  reflection  coefficient  is 
real.  Shonld  the  reflection  coefficient  be  complex  ns  would  be  the  ease  for  total  re¬ 
flection,  then  the  reflected  field  consist.;  of  two  parts;  the  first  is  proportional  to/(f), 
and  the  time  variation  of  the  second  is  given  by  ( Friedlandcr,  194S) 

^rrcjiP£  rr-1/”  ■  <‘»> 


where  X  is  a  parameter  related  to  the  complex  reflection  coefficient.  For  total  re- 
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flection,  the  eorrection  (10)  iepi'c--eiits  the  ■signal  distortion  v  hen  the  simple  geo¬ 
metric  ray  theorj  (nils  te  apply.  However,  this  contrilmtiou  while  imjiortant  for 
theoretical  reasons  turns  out  to  he  a  comparatively  minor  pert urhat ion.  (Arons  ct 
at,  19'*0;  Cron  rf  nl,  I!MW»).  in  nut  analyse.,  we  shall  lind  that  relatively  few  events 
corres|K<iul  to  total  reflection  in  addition,  we  are  le*s  interested  in  calculating  one 
seismogram  exactly  than  in  discerning  the  o|ierations  that  carry  the  response  of 
one  station  into  another,  so  that  small  deviations  common  to  all  traces  will  not  con¬ 
cern  us.  Thus,  although  we  could  easily  correct  the  response  to  account  for  total  re¬ 
flection  by  including  the  contributions  (10)  we  have  not  done  so  for  purjwscs  of 
simplicity.  If  an  event  corrcsjxmds  to  total  reflection,  then  wo  force  the  reflection 
coefficient  to  be  ±1  dcj>ending  upon  its  value  at  the  critical  angle. 


Computation.  Sup|>osc  the  incident  P  wave  lias  the  Unit  variation /(f)  at  soijie 
fixed  station  on  the  crust’s  surface.  It  is  importuut  to  note  that  within  the  theory  of 
geometric  acoustics  the  time  variation  of  /(/)  ami  its  derivative  f  (J)  will  be  pro* 
served  for  all  the  multiple  reflect  cal  rays  except  for  a  translation. 

In  the  acquel,  we  shall  choose /(f)  as 


t  £  0 
t  <  0. 


(ID 


The  virtue  of  this  expression  is  that  with  a  pro|>or  choice  of  the  parameters  o,  0,  y, 
just  about  any  realistic  pressure  variation  of  the  incident  wave  can  be  modeled. 
Typically,  /(f)  =  0  for  f  negative,  thereafter /(f)  rises  rapidly  to  some  maximum 
value,  and  then  decays  to  zero  as  f  becomes  moic  jrositivc  (Figure  o).  It  all  our 
seismograms  wc  have  chosen  tlicvalucso  =.1,0  =  1/4, and y  =  2.  Wc shall  not  need 
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/(<)  explicitly,  but  rather  its  derivative  /(f)  «*  df/dt  which  is  stored  in  the  com¬ 
puter  in  digitised  form  (Figure  G). 

From  (7)  we  calculate  that  the  ground  motion  at  n  station  A  corre8|»nding  to 
any  event  e»\  •••#§')  At  station  .1  is  pro|>ortionn!  to 

u J(t)  m  j" d*  **  )  — ■  t*y/\  —  (e»'),l /(t  ~  r.  ) 


si  to  - - «  w  «  tr  as  as 

Fio.  ft.  A  tree  distrain  and  seismograms  for  an  angle  of  incidence  of  +30*.  In  this  case,  a 
partial  tree  diagram  has  been  drawn  by  impressing  any  ray  whose  amplitude  faila  to  exceed 
.005  times  that  of  the  largest  event  amplitude  in  storage.  The  pressure  variation  of  the  inei* 
dental  wave  front  ia  indicated  on  the  tree  diagram.  This  insert  is  normalised  to  a  scale  appro¬ 
priate  for  the  central  station. 


if  the  incident  wave  is  a  P  wave  (<i  a  positive  integer),  otherwise  we  have 

He'd)  -  [vi  -"(*)’(*/  -  eX  + *  e,v]/(f  -  t,’) 
r.'(f)  -  -  [«.*(«/  +  ti)  +  5  e«Vl  -  (•.')*]/«  “  u) 


(13) 


if  an  SV-wave  is  incident  («i  a  negative  integer).  We  note  that  the  site  of  any  station 
affects  the  response  only  through  the  term/ (<  —  u),  that  is,  a  variable  translation. 
The  complete  seismogram  is  then  found  by  a  superposition  of  these  basic  events. 
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Thu*,  the  left-right  motion  nt  *tntiou  A  i*  given  by  ^  w/(0» 

•mti 

niid  the  U|mIo\vi»  motion  if  given  by  Y  r/OL 

Ml 


(14) 

(15) 


Likewise,  if  we  substitute  6  or  c  for  Motion  ft  or  C’s  roordinntc  we  can  compute  their 
seismograms  too. 


i 


Fig.  7.  A  tree  diagram  and  seismograms  lor  an  angle  of  incidence  of  0*.  In  this  case,  a  partial 
tree  diagram  has  been  drawn  by  mpresaing  any  ray  whose  amplitude  fails  to  exceed  .005  times 
that  of  the  largest  event  amplitude  in  storage.  The  pressure  variation  of  the  incidental  wave 
front  is  indicated  on  the  tree  diagram.  This  insert  is  normalised  to  a  scale  appropriate  for  the 
rent ra'  station. 


C.  Sample  Seismograms.  Wc  have  repeatedly  emphasized  that  the  geometry  of 
the  wedge  problem  is  scale  invariant.  However,  the  time  variation  of  the  pulse  intro¬ 
duces  a  significant  parameter  into  the  analysis.  Consider  any  time  variation  of  the 
incident  wave.  Clearly,  sufficiently  near  the  vertex  it  will  scent  like  a  pulse  of  long 
duration,  and  far  away  front  the  vertex  it  will  have  the  characteristics  of  a  shnrp 
pulse.  Let  fo  be  the  time  required  for  the  pulse  to  build  up  from  its  initial  value  to  its 
maximum.  For  /(f)  given  by  equation  (11)  and  a  «  1, 0  —  J,  nttd  r  -  2,  this  will 
set  ft  to  be  1.2  time  units.  If  wc  multiply  f«  by  the  P  velocity  in  the  mantle  wc  obtain 
a  figure  of  3.6  distance  units.  Wc  can  divide  this  distance  unit  by  the  distance  of  the 
station  from  the  edge  of  the  vertex.  This  dimensionless  number  which  we  shall  write 
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ts  *  is  u  unit  which  indexes  the  scisniogi'nms.  In  detail,  we  have 

/«-iy  nuuitle 
station  distance’ 

There  arc  two  extreme  cases  when  iiitersymlwl  distortion  will  not  occur:  if  / (<) 
is  ei  tiler's  delta  function  (*  «■  0)  or  a  unit  step  Junction  (*  -  »).  In  these  Situations 
the  duration  of  the  event  is  cither  infinitesimal  or  unbounded. 

In  Figures  fi  to  S  wo  present  some  sample  seismograms  for  a  crustal  wedge  of  10*. 
The  upper  three  curves  are  the  computations  of  the  vertical  motion  that  three 
stations  located  at  A,  B,  and  C  would  rt<  ord.  Typically,  the  radial  distances  of  A 
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Fio.  8.  A  tree  diagram  and  seismograms  for  an  angle  of  incidence  of  780*.  I*-,  this  caee,  a 
partial  tree  diagram  has  been  drawn  by  suppressing  any  ray  whose  amplitude  fail*  to  exceed 
.005  times  lhai  of  tb-  largest  event  amplitude  in  storage.  The  pressure  variation  of  the  inci¬ 
dental  wave  front  ia  indicated  on  the  tree  diagram.  Thia  insert  ia  nonnaliied  to  a  acale  appro¬ 
priate  for  the  central  station. 

and  C  have  been  cnoson  to  be  fuel  on  of  .9  and  l.l  respectively  of  that  of  station  B 
with  a  radial  distance  of  25  uni's.  Likewise,  the  x  values  arc  1.00,  1.44,  and  1.31  re¬ 
spectively.  The  lower  trace  is  a  su|>crpcsitiou  of  each  of  the  preceding  three  trans¬ 
lated  to  a  common  time  origin.  It  will  be  noted  that  except  for  the  first  section  of  ihc 
records,  the  three  traces  separate  and  become  distinct.  As  mentioned  previously  one 
reason  for  this  “rubber  sheet  effect”  tsthceffcctivc  magnification  in  the  timescale  for 
stations  located  at  increasing  distances  from  the  wedge  vertex. 

It  will  be  noticed  that  this  is  but  one  distortion:  another  perturbation  is  the  effect 
of  what  might  be  called  intersymbol  interference.  The  seismogram  is  constructed  by 
the  repeated  superposition  of  a  given  wave  shape /(f)  with  various  translations  of 
initial  time  origin.  Since  /(<)  is  not  linear,  the  shape  of  the  combination  signal  of 
two  contributions  will  vary  as  the  relative  spacing  ’-c  tween  the  events  changes.  Ow¬ 
ing  to  the  rapid  decay  of  field  amplitude  as  the  rays  are  multiply  reflected  from  the 
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elliptic  interfaces,  the  effect  of  inlorsyinliol  distort  ion  is  m  n-r  dilliciill  io  set*  for  rela¬ 
tively  few  events  have  significant  amplitude.  l  or  this  mu-sm,  we  have  produced  nil 
artificial  seismogram  shown  in  1'igure  1).  In  this  ilhi>lrnlio;i,  the  motion  was  com- 
pu'ed  by  assuming  that  the  rays  have  rattled  between  the  walls  of  the  wedge  with  an 
artificial  reflection/conversion  eoclhricnt  of  .9  nt  each  bounce.  In  this  fashion,  the 


Flo.  9.  An  artificial  seismograt  comp  iled  by  selling  alt  reflect  ion/cmiversion  coefficients 
to  .9.  The  angle  of  incidence-  is  -30’  Note  ihc  shape  distortion  owing  to  inlcrsynihnl  inter¬ 
ference. 


time  decay  on  the  theoretical  seismogram  is  ntnrh  less  rapid  and  many  more  events 
contribute  to  the  plot  so  that  the  effect  of  intcrsymbol  distortion  is  emphasized. 

For  the  relatively  low  order  rrnstal  model  of  a  uniform  wedge  overlying  an  ideal 
clastic  haif-sim<-e,*  “rubber  sheet  equalisation"  would  be  a  rather  effective  tech¬ 
nique  for  adjusting  the  response  of  one  station  relative  to  another.  However,  in  some 
extensions  of  this  present  work  we  have  considered  more  realistic  custal  confignra- 
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timm,  c.g.,  stacks  of  skew  imiili-iaycm  ovcilyihg  u  half  space.  l  or  tlie.-s*  models,  the 
greater  con\i>lesily  of  the  mishit  *1  nature  is  i  rflntnl  by  haring  cur  res  finn 'ling  teitmu- 
grams  irhirh  arc  much  richer  in  ihlml.  In  these  cases,  “nihl>er  sheet  equalization”  is 
ineffective  mill  more  sophisticated  techniques  me  reqnitisl. 

CONCLUSIONS 

One  a’tn  of  the  present  paper  Inis  liecn  to  show  that  theoretied  seismograms  can 
easily  lie  roust rnetnl  for  cerium  problems  ilirertly  hi  the  time  domain  when  the 
corresponding  situation  would  lie  entirely  intractable  in  the  Imnuonie  domain. 
Imlecil  the  problem  of  diffraetiou  in  wedge-shaped  clast ie  regions  is  one  that  1ms  yet 
to  yield  to  the  mathematician's  touch.  In  the  present  rase,  the  analysis  is  such  that 
seismograms  can  lie  computed  toi  angular  structures  as  easily  as  for  plane  parallel 
geometries.  Furthermore,  the  actual  vibrations  are  easily  and  directly  computed 
without  the  need  to  sum  and  transform  slowly  converging  residue  series. 

It  is  also  important  to  note  that  simple  superposition  of  the  resiioiise  cl  caeh 
clement  in  an  array  is  an  inefficient  practice  and  useful  only  for  the  first  few  seconds 
of  the  initial  motions.  Even  in  the  complete  absence  of  noise,  there  will  not  be 
perfect  coherence  between  the  elements  of  an  array  when  recording  the  same  source. 
However,  a  "rubber  sheet"  equalization  of  the  seismic  traces  will  improve  the  signal 
enhancement  if  the  array  is  located  c  .  ~  dipping  crustal  structure.  On  tlm  other 
hand,  this  rubber  sheet  equalization  will  not  compensate  for  the  intersymbol 
distortion.  However,  neither  of  these  practices  are  necessary,  for  if  theoretical 
seismogram  can  be  constructed  for  a  realistic  earth  model  by  linear  supemosition, 
then  it  is  relatively  easy  to  construct  a  “matched  filter"  which  run  greatly  enhance 
Ihe  signal/noise  ratio.  In  subsequent  pajiers  wc  plan  to  illustrate  the  computation 
of  theoretical  se  is. nog  rains  for  skew  multilayers,  and  wc  shall  discuss  the  construc¬ 
tion  of  matched  filters  1  ssed  upon  this  analysis. 
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(p"rl  r  a  series  or;  crustal  equalization  of  seismic  arrays) 


by  Julius  Kane  and  Dr  or  Gan  oh 


abstract 

Most  seismic  theory  is  confined  to  crustal  structures  that  can  be 
w  is t. -acted  as  plane  parallel  layers.  For  such  configurations,  the  response 
of  cnn  element  of  a  seismic  array  will  be  similar  to  any  other  except  for 
a  tine  delay.  Therefore  signal  enhancement  for  such  an  array  can  be  effect¬ 
ed  quite  ea  lly  by  suitable  time  delays  of  the  individual  records  followed 
by  superposition  to  reduce  the  effects  of  random  noise. 

On  the  other  hand,  if  the  crustal  structure  i-  anything  but  a  plane 
parallel  configuration,  the  signals  received  by  the  individual  elements  will 
not  be  translated  copies  of  one  another,-  but  will  include  distortion 
effects  characteristic  of  the  local  geometry.  As  a  result,  the  signals 
received  by  seismic  arrays  located  on  realistic  crustal  configurations 
need  to  be  equalized  to  some  standard  reference  if  optimum  signal  proces¬ 
sing  is  to  be  achieved. 

In  this  paper  we  describe  a  procedure  for  the  calculation  of 
theoretical  seismogram?  for  the  tel ose-irznic  response  of  an  array  of 
stations  located  above  a  non-uniform  crust  composed  of  skew- layers,  xn 
terms  of  our  mathematical  model  ve  demonstrate  the  distortion  effects  of 
such  a  structure  and  discuss  equalization  techniques  that  will  permit  a 
superio:  recovery  of  the  desired  signal. 


*  Presented  at  the  IUGG  symposium  on  geophysics  ar.d  computers  at 
Cambridge  University,  Cambridge,  England,  June  i960. 
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A  CCKP7TATI0NAI  APPROACH  TO  TfE  INVERSE  STRAIN  PROBLEM* 

1 

by  Julius  Kane  and  A jit  Kul 
Space  Science  Center 
University  of  California 
Los  Angeles,  California  5002^ 


abstract 


In  this  paper  we  derive  algorithms  for  the  solution  of  the 
following  problem:  An  elastic  half  space  is  deformed  from  its 
rest  position  by  some  internal  force  distribution.  The  boundary 
displacement  is  assumed  to  be  known.  The  task  posed  is  to  find 
a  suitable  force  distribution  in  the  interior  of  the  solid  which 
can  account  for  the  strain  of  the  boundary. 

We  solve  the  problem  in  the  following  fashion:  We  introduce 
a  pair  of  functions  which  have  che  property  that  they  are  solutions 
of  the  relevant  partial  differential  equations  and  can  be  used  to 
generate  a  double  sequence  of  functions  which  can  interpolate  the 
vector  data  on  the  boundary.  Thus,  we  show  that  if  an  arbitrary 
point  is  chosen  as  the  epicenter  of  the  force  distributor ,  then  a 
source  field  centered  at  that  point  can  account  for  any  surface 
displacement,  to  any  desired  degree  of  precision.  The  ambiguity 
stems,  in  part  from  the  lack  of  uniqueness  in  the  geophysical  problem, 
but  more  frora  the  instability  of  the  solution  with  respect  to  the 
boundary  data:  That  is,  a  small  change  in  the  surface  displacement 
can  imply  a  large  change  in  the  source  distribution. 

^Presented  at  the  IUGG  Symposium,  Cambridge,  England,  June  19 66. 
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abstract 

Consideration  is  given  to  a  mathematical  analysis  of  a  special 
case  of  seismic  vave  phenomena.  In  particular,  the  mined  boundary 
value  problem  of  Love  wave  propagation  in  a  solid  layer  over  a  solid 
half-space  is  investigated  where  the  layer  undergoes  an  abrupt  change 
ir.  thickness.  Both  the  layer  and  half-space  arc  considered  to  be 
homogeneous  clastic  media.  Theoretical  background  fo  the  physics  of  . 
the  problem  is  provided  by  presenting  the  fundamentals  of  clastic 
wave  propagation  with  specialisation  to  the  Love  wave  case  ar.c.  by  a 
statement  of  the  physical  nature  ar.d  mathematical  form  of  applicable 
boardary  conditions.  Interest  is  focused  or.  the  amplitudes  of  the 
transmitted  and  reflected  Love  waves  relative  to  the  magnitude  of  the 
excitation,  i.e.,  the  transmission  and  reflection  coefficients,  for  a 
range  of  both  the  "strength"  of  the  discontinuity  (magnitude  of  change 


in  thirlviosn  relative  to  layer  thickncuc)  and  the  layer  thickness. 

The  analysis  employs  scattered  fields  in  the  form  of  integrals  in 
the  complex  plane;  boundary  conditions  are  applied  to  the  total  fields 
lb  satisfy  the  resulting  equations  for  the  boundary  conditions,  the 
nature  of  some  of  the  unknown  coefficients  in  the  scattered  field 
integrals  is  postulated  and  a  function- theoretic  argument  is  employed 
to  determine  these  coefficients.  The  transmission  end  reflection 
coefficients  are  then  extracted  by  a  standard  appeal  to  the  calculus 
of  residues  and  the  energy  contained  in  the  diffracted  field  is 
evaluated.  Both  displacement  amplitude  and  energy  coefficients  are 
displayed  graphically  as  a  function  of  the  parameters  h  and  H.  The 
results  presented  show  that  the  reflected  energy  constitudes  less 
than  one  percent  of  the  incident  energy  for  all  crust  thicknesses 
considered  here.  The  amplitude  transmission  coefficient  relating 
the  relative  magnitude  of  the  displacements  of  the  transmitted  and 
incident  waves  is  shown  to  take  on  values  greater  than  unity  for  low 
frequency  and  an  energy  analysis  shows  that  this  behavior  does  r.ot 
violate  the  principle  ?f  energy  conservation.  A  comparison  reveals 
that  the  results  shown  here  are  in  close  agreement  with  those  given 
by  Knopoff  and  Hudson  (reference  8).  However,  a  discrepancy  does 
exist  for  intermediate  values  of  crust  thickness  and  is  discussed  in 


some  detail. 
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